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ABSTRACT 


In  this  investigation  the  general  solution  is  derived  for  the 
problem  of  the  optimum  linear  estimation  of  a  sampled  stochastic 
process,  when  the  transition  and  output  matrices  of  the  model  of  the 
process  are  random  parameters  that  are  independent  from  one  sample  point 
to  the  next  with  known  mean  and  covariance.  The  resulting  estimate  is 
optimum  in  the  sense  that  it  minimizes  the  trace  of  the  covariance  matrix 
of  the  error  (a  generalized  mean-squared-error  criterion). 

The  notation  used  in  the  following  discussion  is  based  on  the  state- 
transition  approach  to  linear  estimation  developed  by  Kalman,  In  this 
approach  the  stochastic  process  is  represented  as  the  output  of  a 
linear  (possibly  time-varying)  dynamic  system  with  an  independent  random 
input. 

For  current  estimation  and  prediction  of  the  state  vector,  the 
optimum  estimate  is  implemented  by  a  linear  dynamic  filter  with  a 
matrix-valued  gain  the  only  undetermined  coefficient.  This  matrix¬ 
valued  gain,  as  well  as  the  covariance  matrix  of  the  error  in  the 
optimum  estimate,  is  determined  iteratively  for  each  sample  point  from 
a  nonlinear  difference  equation  involving  the  covariance  of  the  error 
at  the  previous  sample  point. 

The  configuration  of  the  solution  for  linear  interpolation  with 
delay  is  a  linear  dynamic  filter  similar  to  the  one  used  for  prediction. 
For  each  sample  period  Che  estimate  is  delayed,  an  additional  weighting 
matrix  and  delay  element  must  be  added  to  the  filter. 

All  of  these  results  are  derived  from  the  sampled  version  of  the 
Wiener-Hopf  equation,  and  they  apply  without  modification  to  stationary 
and  nonstationary  statistics  and  to  growing-memory  and  infinite-memory 
filters . 
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I  INTRODUCTION 


A.  STATEMENT  OF  THE  PROBLEM 

This  investigation  concerns  the  linear  estimation  of  a  sampled 
stochastic  process  when  certain  parameters  of  the  process  must  be  treated 
as  random  variables  that  are  independent  from  one  sample  point  to  the 
next.  These  random  parameters  may  be  due  to  multiplicative  noise  that 
is  corrupting  the  observed  samples,  to  random  variations  in  the  sample 
period,  or  to  other  uncertainties  in  the  a  priori  knowledge  of  the  pro¬ 
cess.  The  stochastic  process  is  represented  (for  average  statistical 
properties  up  to  second  order)  as  the  output  of  a  linear  dynamic  system 
excited  by  independent  gaussian  processes.  This  model  of  a  stochastic 
process  is  very  general,  and  in  particular  it  includes  the  important 
special  case  of  stationary  statistics  and  rational  power  spectra  as  well 
as  a  large  class  of  nonstationary  processes.  Both  the  stochastic  process 
and  the  random  parameters  may  be  stationary  or  nonstationary,  and  the 
linear  estimation  includes  the  case  in  which  the  number  of  observed 
samples  is  growing . 

By  way  of  example,  consider  the  following  problem.  A  satellite  is 
telemetering  data  to  a  distant  ground  station  The  original  data  are 
continuous,  but  they  must  be  sampled  before  they  are  transmitted.  Noise 
in  the  electronics  or  random  fading  in  the  transmission  characteristics 
of  the  atmosphere  can  introduce  multiplicative  noise.  The  time  between 
successive  sample  points  may  vary  in  a  random  manner  because  of  jitter 
or  missed  samples  in  periodic  sampling  resulting  from  imperfections  in 
the  equipment,  jamming,  or  natural  interference.  On  the  other  hand,  the 
data  may  be  transmitted  at  random  intervals  intentionally  because  of  the 
random  character  of  the  quantity  being  measured  or  in  order  to  counter¬ 
act  jamming.  In  other  words,  the  ground  station  may  be  operating  on 
randomly  sampled  data  with  multiplicative  noise.  From  these  observed 
samples  it  is  desired  tc  obtain  a  continuous  estimate  of  the  original 
data  or  to  predict  the  value  of  the  data  at  some  future  time. 

The  output  of  a  communication  system  such  as  the  one  described  here 
can  be  represented  as  a  stochastic  process,  A  realistic  model  of  the 
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stochastic  process  must  include  all  the  sources  of  random  variation.  In 
this  investigation  the  model  of  the  stochastic  process  is  based  on  a 
combination  of  the  Bode-Shannon  representation  of  a  random  process  and 
the  "state-transition"  method  of  analysis  of  dynamic  systems  introduced 
by  Kalman  TRef„  1],  The  output  matrix  and  the  transition  matrix  of  the 
sampled  model  are  matrix-valued  random  parameters  because  of  multiplica¬ 
tive  noise  and  random  variations  in  the  sample  period.  The  random  para¬ 
meters  have  a  known  probability  distribution  (not  necessarily  stationary) 
independent  from  one  sample  period  to  the  next. 

The  optimum  estimate  of  the  state  variables  is  the  linear  estimate 
which  minimizes  the  trace  of  the  covariance  matrix  of  the  error.  This 
is  a  natural  extension  of  the  linear  least-squares  estimate  discussed 
extensively  in  the  literature  [Ref,  2]  If  the  stochastic  process  is 
actually  generated  by  gaussian  random  excitations  of  a  linear  dynamic 
system,  and  if  the  probabllit5r  distribution  of  the  random  parameters  is 
unimodal  and  symmetric  about  the  sxpected  value  of  the  parameters,  then 
the  optimum  linear  estimate  gives  the  conditional  expectation  of  the 
desired  state  variables.  Sherman  [Ref.  3]  has  shown  that  this  condi¬ 
tional  expectation  minimizes  the  expected  value  of  a  large  class  of  loss 
functions,  Gunckel  [Ref.  4]  has  proved  that,  when  the  state  variable  is 
not  known  exactly,  its  conditional  expectation  can  be  used  in  the  solu¬ 
tion  to  the  general  problem  of  control  with  a  quadratic  loss  function. 
Therefore,  one  is  led  to  the  conclusion  that  the  conditional  expectation 
is  the  best  estimate  in  the  general  control  problem  as  well  as  in  the 
estimation  problem.  On  the  other  hand,  the  optimum  linear  estimate 
requires  only  the  average  statistics  of  the  process  up  to  second  order, 
and  if  only  the  mean  and  the  covariance  of  the  process  are  known,  then 
it  is  the  best  estimate  that  can  be  made  with  this  information. 

In  addition  to  estimating  the  current  value  of  the  state  variables, 
it  may  be  necessary  to  predict  the  value  of  the  state  variables  at  some 
time  in  the  future,  or  it  may  be  advantageous  to  interpolate  some  past 
values  of  the  state  vector  from  more  recently  observed  random  variables. 
The  interpolation  with  delay  should  reduce  the  trace  of  the  covariance 
of  the  error  because  more  random  variables  have  been  observed  during  the 
delay. 
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B.  PREVIOUS  WORK 


The  early  work  of  Wiener  [Ref.  5]  showed  that,  for  a  continuous 
stochastic  process,  the  problem  of  linear  estimation  leads  to  the 
Wiener-Hopf  integral  equation.  For  the  practically  important  case  of 
stationary  statistics  and  rational  power  spectra,  he  demonstrated  that 
the  solution  to  the  integral  equation  could  be  obtained  by  spectral 
factorization.  Under  the  same  conditions  of  stationary  statistics  and 
rational  power  spectra,  Franklin  [Ref,  6]  solved  the  problem  of  linear 
estimation  using  periodically  sampled  data,  and  Amara  [Ref.  7]  general¬ 
ized  this  work  to  include  multivariable  systems. 

Concurrently,  there  has  been  considerable  interest  in  the  literature 
on  the  analysis  and  stability  of  systems  with  random  parameters  [Ref.  8, 
9],  although  not  much  work  has  been  concerned  with  the  design  of "filters 
for  these  systems.  Kalman  [Ref-  10]  considered  the  optimum  control  of 
a  linear  system  that  is  randomly  sampled  using  a  quadratic  error  cri¬ 
terion  and  a  step  input,  Gunckel  [Ref,  U]  has  extended  more  recent  work 
of  Kalman  [Refs.  11,  12]  to  provide  a  general  solution  to  the  problem 
of  the  control  of  a  linear  system  with  random  parameters .  In  particular, 
he  shows  that,  if  it  is  desired  to  minimize  the  expected  value  of  a 
quadratic  loss  function,  the  conditional  expectation  of  the  state 
variables  should  be  used  in  the  optimum  control  procedure.  Gunckel 's 
work  separates  effectively  the  problem  of  estimation  of  the  state 
variables  from  the  problem  of  control  of  the  state  variables  His 
results  are  especially  important  because  under  the  conditions  discussed 
in  the  previous  section,  the  optimum  linear  estimate  derived  in  the 
present  investigation  is  the  conditional  expectation  of  the  state 
variable  and  therefore  can  be  used  in  the  optimum  control  procedure. 

For  a  randomly  sampled,  stationary  stochastic  process,  Bergen 
[Ref.  13]  has  determined  the  spectral  density  from  a  convolution  inte¬ 
gral  that  can  be  evaluated  in  some  cases  by  the  method  of  residues.  A 
synthesis  procedure  to  determine  the  best  linear  time -invariant  continu¬ 
ous  filter  for  these  cases  is  based  upon  the  standard  Wiener  factoriza¬ 
tion  of  the  sampled  power  spectrum  [Ref.  14]. 

Buetler  [Refs.  15,  l6]  has  generalized  the  Wiener  theory  to  include 
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stochastic  processes  with  random  parameters  for  continuous  stationary- 
systems.  He  obtained  practical  solutions  for  two  cases — the  optimum 
linear  filter  when  the  prediction  time  or  lag  time  is  a  random  parameter, 
and  the  optimum  linear  estimate  when  the  system  gain  is  a  random  para¬ 
meter  (multiplicative  noise).  Shaw  [Ref,  17]  has  considered  the  problem 
of  dual-mode  filtering  for  a  continuous  stationary  stochastic  process 
when  the  instantaneous  model  of  the  process  varies  randomly  between  two 
possible  models. 

Kalman  [Ref.  1]  has  formulated  the  whole  problem  of  linear  estima¬ 
tion  from  sampled  data  in  matrix  notation  in  terms  of  state -transition. 
The  problem  is  approached  from  the  paint  of  view  of  conditional  expecta¬ 
tions  rather  than  from  the  sampled  version  of  the  Wiener-Hopf  equation. 
Kalman  and  Bucy  [Ref.  18]  have  extended  this  formulation  to  the  linear 
estimation  of  a  continuous  stochastic  process  when  white  noise  is  added 
to  the  measurements.  In  this  extension  the  results  are  derived  from  the 
continuous  Wiener-Hopf  integral  equation. 

The  work  in  this  investigation  represents  a  generalization  of  this 
approach  to  the  problem  of  linear  estimation  in  the  presence  of  random 
parameters . 

C.  OUTLINE  OF  NSW  RESULTS 

The  solution  to  the  general  problem  of  the  optimum  linear  estimate 
of  a  sampled  stochastic  process  with  random  parameters  is  derived  in 
this  investigation  Chapter  II  is  a  review  of  the  state-transition 
approach  to  linear  estimation  using  matrix  notation,  with  appropriate 
examples  of  random  parameters  included. 

In  Chapter  III  the  optimum  filter  for  the  current  estimation  and 
prediction  of  the  process  described  in  the  previous  chapter  is  derived 
from  the  sampled  Wiener-Hopf  equation;  the  desired  weighting  coefficients 
and  the  covariance  of  the  error  in  estimation  are  determined  iteratively 
for  each  sample  point  from  the  covariance  of  the  previous  estimate. 

These  results  are  extended  in  Chapter  IV  to  the  problem  of  optimum 
interpolation  with  delay  (at  time  t  the  optimum  estimate  is  desired 
of  the  state  vector  at  time  t  where  t  <  t  ).  Chapter  IV  presents 
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the  first  thorough  investigation  in  the  literature  of  optimum  interpola¬ 
tion  with  delay  for  a  nonstationary  stochastic  process. 

Chapter  V  applies  the  ideas  developed  in  this  investigation  to  the 
estimation  of  a  stationary  stochastic  process  with  a  random  sample 
period,  and  for  a  simple  example  the  optimum  filter  is  compared  with  the 
best  linear  time-invariant  filter. 

For  current  estimation  and  prediction,  the  optimum  estimate  is 
implemented  by  a  linear  dynamic  filter.  The  only  undetermined  coeffi¬ 
cient  of  the  filter  is  the  matrix-valued  gain,  which  is  determined 
iteratively  for  each  sample  point.  When  the  statistics  of  the  process 
are  stationary,  the  matrix-valued  gain  approaches  a  steady-state  value 
as  the  number  of  sample  points  approaches  infinity 

The  configuration  of  the  optimum  solution  for  linear  interpolation 
with  delay  is  shown  to  be  a  linear  dynamic  filter  similar  to  the  one 
used  for  prediction;  but,  for  each  sample  period  the  estimate  is  delayed, 
an  additional  weighting  matrix  must  be  determined. 

In  order  to  relate  this  investigation  to  the  conventional  approach 
to  linear  estimation,  all  of  these  results  are  derived  from  the  sampled 
version  of  the  Wiener-Hopf  equation 
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II.  DESCRIPTION  OF  THE  SYSTEM 

A.  INTRODUCTION 

This  chapter  will  provide  an  introduction  to  the  state-transition 
approach  to  linear  estimation.  The  state  of  the  system  is  some  quanti¬ 
tative  information  (such  as  a  set  of  numbers)  which  is  the  minimum  amount 
of  data  necessary  to  predict  the  future  behavior  of  the  system.  When 
applied  to  estimation,  the  state  is  the  data  necessary  to  predict  the 
expected  value  of  the  future  behavior  of  the  system.  The  state  transition 
specifies  the  dynamics  of  the  system — how  the  state  at  one  instant  of 
time  is  transferred  to  the  state  at  a  later  instant  of  time. 

This  paper  will  consider  only  those  stochastic  processes  which  can 
be  represented  (for  average  statistical  properties  up  to  second  order) 
as  the  output  of  a  linear  dynamic  system  excited  by  independent  gaussian 
processes.  Therefore,  at  any  instant  of  time,  the  state  of  the  system 
can  be  represented  by  an  r-dimensional  vector,  and  the  state  transition 
is  an  r  x  r  matrix  with  the  properties  enumerated  in  Section  II-C.  This 
representation  is  a  very  general  one,  and  in  particular  it  includes  the 
important  special  case  of  stationary  statistics  and  rational  power 
spectra,  as  well  as  a  large  class  of  nonstationary  processes.  In 
Sections  II-B  and  III-C  this  representation  is  presented  for  two  common 
stationary  processes. 

The  sampled  stochastic  process  is  represented  in  Section  II-D  as  a 
linear-difference  equation  with  random  parameters  and  random  excitation. 
Examples  of  random  parameters  are  discussed  in  Section  II-E  for  the 
output  matrix  (caused  by  multiplicative  noise)  and  the  transition  matrix 
(when  the  sample  period  is  an  independent  random  variable).  It  should 
be  emphasized  that  the  random  parameters  are  not  the  only  source  of 
noise;  both  correlated  and  uncorrelated  noise  can  be  included  in  the 
vector  representing  the  state  of  the  system. 

Finally,  in  Section  II-F  are  discussed  the  reasons  for  restricting 
the  optimum  estimate  to  a  linear  combination  of  the  observed  random 
variables . 

T 

In  the  notation  convention  followed  here,  A  is  the  transpose  of 
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A,  and  E  [A]  is  the  expected  value  of  A.  When  v  is  an  m  X  1 

T 

vector,  then  v  v  is  the  scalar  product  resulting  in  a  scalar,  and 
T 

w  is  the  vector  product  resulting  in  an  m  X  m  matrix.  The  elements 


of  the  A  matrix  are  denoted  by 
by  v±. 


and  the  components  of  the  vector 


B.  MODEL  OF  THE  CONTINUOUS  PROCESS 


The  model  of  the  stochastic  process  is  based  on  the  Bode-Shannon 
representation  of  a  random  process  and  the  "state-transition"  method  of 
analysis  of  dynamic  systems  introduced  by  Kalman  [Ref.  1].  A  linear 
dynamic  system  can  be  described  by  the  following  ordinary  differential 
equation, 


d 

—  x(t)  =  F(t)  x(t)  +  G(t)  v(t) 
dt 

y(t )  =  M(t)  x(t) 


(i) 


where 

x(t)  is  an  r  X  1  vector  that  is  the  state  of  the  system, 

v(t)  is  an  m  X  1  vector  that  is  the  input  uo  the  system, 

y(t)  is  a  p  X  1  vector  that  is  the  output  of  the  system, 

F(t)  is  an  r  X  r  matrix  representing  the  dynamics  of  the  system, 

G(t)  is  an  r  X  m  distribution  matrix  representing  the  constraints 


on  the  input , 

M(t)  is  a  p  X  r  output  matrix. 


The  components  x^  of  the  state  vector  x  are  called  the  state 
variables,  while  the  components  of  the  output  y  are  linear  combinations 
of  these  state  variables.  The  matrix  F  may  be  nonsingular  and  repre¬ 
sents  the  dynamics  of  the  system.  In  stationary  systems  the  matrices 
F,  G,  and  M  are  constants. 

The  matrix  block  diagram  of  the  system  is  presented  in  Fig.  1.  The  • 
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FIG  1,  MODEL  OF  THE  CONTINUOUS  STOCHASTIC  PROCESS. 

thick  lines  indicate  vector-signal  flow,  and  the  transfer  function  l/s 
actually  stands  for  r  integrators  such  that  the  output  of  each  is  a 
state  variable.  The  dynamics  matrix  F(t)  indicates  how  the  outputs 
of  the  integrators  are  fed  back  into  the  inputs  of  the  integrators. 

Thus  f  (t)  is  the  coefficient  with  which  the  output  of  the  jth 
integrator  is  fed  back  to  the  input  of  the  ith  integrator. 

This  linear  dynamic  system  represents  a  stochastic  process  when  the 
input  to  the  system,  v(t),  is  a  random  process.  In  the  model  used  in 
this  investigation,  the  input  v(t)  is  an  m  X  1  vector-valued  random 
process  with  zero  mean  and  with  the  m  X  m  covariance  matrix 

E[v(t)  vT(s)]  =  V(t)  &(t  -  s),  (2) 

where  &  is  the  Dirac  delta  function.  The  defining  property  of  the 
delta  function  is,  for  any  function  V(t)  bounded  and  continuous  at  s, 

8 (t  -  s )  =  0  t^s 
00 

^  V(t)s(t  -  s)dt  =  V(s ) .  (3) 

-00 

This  definition  is  satisfactory  because  only  the  integral  of  the  covar¬ 
iance  in  Eq.  (2)  is  ever  required-  Because  the  input  v(t)  is  an  inde¬ 
pendent  random  variable  with  zero  mean, 
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E[v(t)  xT(t0)]  =  0 


t  >  tQ.  (k) 

When  the  stochastic  process  is  stationary  with  a  rational  power 
spectrum  S  (s),  the  transfer  function  H(s)  of  the  linear  dynamic 
system  representing  the  process  is  obtained  by  factoring  the  power 
spectrum 


S  (s)  =  H(s)  H(-s)  (5) 

where  H(s)  contains  all  the  poles  and  zeros  in  the  left  half  plane. 

For  example,  consider  the  first  order  stationary  Markov  process 
with  the  power  spectrum  and  autocorrelation  function 


S  xx(s) 


2P 


Q2  2 
P  -  S 


$  (r)  =  e"plTI 

°  vv  ' 


XX 


(6) 


The  power  spectrum  is  the  Fourier  transform  of  the  autocorrelation 
function  [Ref.  19]}  therefore,  both  functions  give  the  same  information 
about  the  process.  In  this  example  the  transfer  function  of  the  linear 
system  is 


H(s) 


(ap)1/g 

p  +  s 


(7) 


The  linear  differential  equation  describing  the  process  is 


=  -  Px(t)  +  v ( t ) ,  (8) 

with  the  input  v(t)  an  independent  random  variable  with  zero  mean  and 
covariance 


E[v(t)  v(s)]  =  vlx  6 (t  -  s)  =  2p  b(t  -  s).  (9) 

The  model  of  the  stationary  Markov  process  is  presented  in  Fig,  2. 
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v(t) 


FIG.  2.  MODEL  OF  A  STATIONARY  MARKOV  PROCESS. 


The  general  solution  to  the  ordinary  differential  equation,  Eq.  (l). 


is 


x(t)  =  $(t,t0)x(tQ)  +  I  f (t,T)G(r)v(T)dr  t  >  tn,  (10) 


-  O’ 


where  $ (t,tQ)  is  an  r  X  r  matrix  called  the  transition  matrix  of 
the  system.  From  Eqs.  (4)  and  (10),  it  is  easily  seen  that 

E[x(t)xT(tQ)]  =  i(t,tQ)  E[x(t0)xT(t0)]  t  >  tQ,  (11) 

since  x(t^)  is  independent  of  v(t)  for  t  >  t^. 


C.  TRANSITION  MATRIX 


Tue  pertinent  information  about  the  transition  matrix  is  presented 
in  this  section.  For  a  more  complete  treatment  see  Coddington  and 
Levinson  [Ref.  20].  The  transition  matrix  is  a  nonsingular  matrix 
satisfying  the  differential  equation 

d 

—  I  =  F(t) 5  ,  (12) 

dt 


made  unique  by  the  requirement  that 
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(13) 


*  <vV  ■ l’ 


where  I  is  the  identity  matrix. 

Two  properties  of  the  transition  matrix  are 


1  (^3^2^  - 


2*  1 


l’u2' 


(14) 


If  the  matrices  F(t)  and  J  F(r)dT  commute,  then  the  transition 

*0 

matrix  can  he  written 


exp 


F(T)di 


where  the  definition  of 


exp(B) 


B 
k : 


(15) 


(16) 


For  a  discussion  of  the  properties  of  commuting  matrices  see  Gantmacher 
[Ref.  21].  In  particular,  the  two  matrices  commute  if  F(t)  is  diagonal 
or  if  F  is  a  constant.  When  the  F  matrix  is  constant,  the  transition 
matrix  $  is  stationary,  and  it  can  be  written 

5  (tQ  +  T,tQ)  =  $  (t)  =  exp(F-r)  =  eFr,  (17) 

where  exp(Fr)  is  called  the  exponential  of  the  matrix  Ft. 

When  the  characteristic  roots  of  F  are  distinct,  then  the  matrix 
F  is  similar  to  a  diagonal  matrix  A  with  diagonal  . . . ,A  ,  so 

that 

F  =  DAD-1,  (l8) 


11 


where  D  is  a  nonsingular  r  X  r  matrix  and  the  characteristic  roots 


may  be  complex.  Therefore,  the  transition  matrix  is 


$  (t)  =  exp(Fr)  =  D  •  exp(A-r)  •  D  t 


(19) 


This  is  readily  seen  to  be  the  case  because 


m  -l.k  k 

_  V'  (DAD  )  t 

exp(DArD_  )  =  I  +  / 


=  D<  I  + 


k=l 


ki 


l 


(At  Y 


.-1 


k=l 


kl 


(20) 


This  method  of  taking  the  exponential  of  a  matrix  is  very  satisfactory 
when  F  is  already  diagonal  or  nearly  so,  but  sometimes  in  actual  prac¬ 
tice  this  diagonalization  may  be  difficult  to  perform,  and  an  alternate 
method  may  be  more  efficient. 

An  alternative  method  for  taking  the  exponential  of  a  matrix, 
given  by  Friedman  [Ref.  22],  is  based  on  the  following  theorem  in  his 
chapter  on  Spectral  Theory  of  Operators: 

"If  A  is  a  matrix  whose  eigenvalues,  arranged  in  order 
of  increasing  absolute  value  are  A^, Ag, . . . ,An  and  if  g(A) 
is  an  analytic  function  of  A  in  a  circle  around  the  origin 
with  radius  greater  than  | An  | ,  then  g(A)  equals  r(A), 
the  polynomial  of  degree  n  -  1  for  which 

g(\)  =  r(\)  k  =  1,2,..., n."  (21) 

In  particular,  this  means  that  if  A  is  an  n  X  n  matrix  with  distinct 

characteristic  roots  \  ,  then 

1  d  n 

n 

eAt  =  aQI  +y  ajAt)1,  (22) 

i=l 

where  the  a  are  evaluated  by  the  set  of  equations 


12  - 


(23) 


e  —  0:^  +  ^  ^(^t)  k  -  lj 2,  .  ..,n. 

i=l 

When 


so  that  is  a  characteristic  root  of  multiplicity  v  for  A 

Eq.  (21)  is  modified  as  follows: 


gJ(\) 

-  Av 

J  =  0,1,2,.. 

g(\) 

■  r<v 

k  =  v+1,  v+2 

Thus  Eq.  (23)  is  modified  to 


i  =  0,1,2, . . „,v-i 


k  =  v+1, v+2, . . . ,n 


i=l 


To  illustrate  this  method,  the  exponential  will  "be  taken  of 
arbitrary  2X2  matrix  At  where 

A  = 


The  characteristic  roots  are  determined  from 


‘11 


21 


‘12 


22 


|A  -  Al|  =  0 

(&11  '  *)(a22  ‘  -  ai2S21  =  0 


,  then 


(25) 


(26) 

the 


(27) 
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(28) 


X  =  i(alx  +  a22)  +  i  >/(an  -  a^)  +  Ua^ 


The  two  characteristic  roots  A^  and  X2  are  distinct  if 


(ail  *  a22}  +  4&12a21  *  0 


(29) 


For  the  2X2  matrix  Eq.  (23)  is 


v 

e  =  a0  +  a^t 

v 

e  =  ao  +  V2t} 


(30) 


with  the  solution 


ao  = 


At  A  t 

xie  -  V 
\'X2 


Ajt  Agt 

e  -  e 


V  = 


\  -  X2 


(31) 


Thus,  when  the  characteristic  roots  of  the  2X2  matrix  A  are 
distinct,  the  exponential  of  the  matrix  At  is 


A2t  At 

A  e  -  A  e 

.At 

"l 

o' 

A,  t  At 

e  1  -  e  2 

ail 

S12 

0 

1 

al  -  a2 

_S21 

S22_ 

(32) 


If  the  characteristic  roots  are  not  distinct,  then  from  Eq.  (26),  the 
exponential  of  the  matrix  At  is 


At  V 

e  =  e 


ail  +  1  -  V  ai2 
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a22  +  1  -  V 


(33) 


As  an  example,  the  exponential  in  Eq.  (32)  will  be  used  to  calculate 
the  transition  matrix  of  a  stationary  stochastic  process  with  two  state 
variables.  Consider  the  process  with  the  following  power  spectrum  and 
autocorrelation 


v  (P  +  s)(p  -  s)  L 

t>x  X  (s)  =  - P - P - p - p"- 

f(e  +  s)2  +  72][(p  -  s)2  +  7] 


Sx^Ct)  =  e 


•Mr|  cos  yT  _  gz^azlij 
V  2p  +  7 


L  = 


+  72) 

2p2  +  72 


(3M 


The  transfer  function  of  the  dynamic  system  is 


H(s) 


=  (P  +  *0  L1^2 
(P  -  s)2  +  y2 


(35) 


The  model  of  this  stochastic  process  is  presented  in  Fig.  3,  where  the 
two  state  variables  are  x^(t),  the  observed  random  variable  and  the 
output  of  one  integrator,  and  x2(t) ,  the  output  of  the  other  integrator. 
The  two  linear  differential  equations  describing  the  system  are 


v ,  ( t) 


FIG.  3.  MODEL  OF  A  PROCESS  WITH  AN  EXPONENTIAL  COSINE  AUTOCORRELATION 
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and 


dx  (t) 

-  =  -  Pxx(t)  +  x2(t)  +  v1(t) 

dt  “ 


<ix2(t)  0 

- -  -  7^00  -  ex2(t). 

dt 


The  input  v^t)  has  zero  mean  and  covariance 

E[v1(t)v1(s)]  =  vxl  6 (t  -  s). 
For  this  system  the  2X2  dynamics  matrix  F  is 

F  = 


and  the  characteristic  roots  are  complex  with 


\  =  -  P  -  *7 
=  "  P  +  J7 

where 

j  is  the  imaginary  unit  number. 

From  Eq.  (32)  the  exponential  of  the  matrix  Ft  is 


*(t)  = 


Ft 

e 


-Pt 


cos  7t  1/7  sin  7t 
-  7  sin  7t  cos  7t 


where  the  sum  of  the  complex  conjugate  exponentials  has  "been 
terms  of  sines  and  cosines 


(36) 


(37) 


(38) 


(39) 


(40) 
written  in 
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sin  yt  =  -=j  j  (e  -  e+^^) 

cos  yt  =  -g  (e-^  +  e+j7t)  (41) 

The  model  of  the  first-order  stationary  Markov  process  in  Fig.  2  and  the 
model  of  the  process  with  exponential  cosine  autocorrelation  in  Fig.  3 
will  he  used  in  the  numerical  examples  calculated  in  the  following 
chapters. 

D.  MODEL  OF  THE  SAMPLED  PROCESS 

In  many  cases,  a  linear  dynamic  system,  such  as  that  shown  in  Fig.  1, 
is  sampled  at  points  in  time  tQ,t  , . . .t  , . . .t  called  the  sample 
points,  where  the  subscript  k  indicates  the  kth  sample.  The  time 
between  successive  samples  is  called  the  sample  period  T  ,  which  is 

.K. 

given  by 


Tk  = 


“k+1 


V 


(42) 


It  is  assumed  that  all  the  switches  in  the  sampling  operation  operate 
synchronously  and  that  the  sampling  operation  can  be  represented  as 
the  result  of  modulating  an  impulse  train  B,p(t)  with  the  output  of 
the  system  y(t)  so  that 


y*(t)  =  y(t)  &T(t),  (43) 

where  the  impulse  train  6,p(t)  is  given  by 

00 

6T(t)  6(t  -  tk).  (44) 

k=0 


When  a  linear  system  with  an  impulse  i'3sponse  h(t)  follows  the  sampling 
operation,  then  the  output  of  the  linear  system  at  time  t  is 
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:(t)  =y  h(t  -  t  )  y(t  ) 


t  <  t  <  t  . 
n  n+1 


(*5) 


k=0 


If  a  clamp  or  zero-order  hold  follows  the  sampling  operation,  then  the 
output  of  the  hold  at  time  t  is  y(t  )  with  t  <  t  <  t 

When  the  state  of  the  system  is  considered  only  at  the  sample 
points,  the  stochastic  process  is  said  to  be  discrete,  and  the  model 
under  consideration  becomes 


x(tn+l)  =  f  (Wtn)  X(tn}  +  u(tn} 


y(tn)  =  M(tn)  x(tn), 


(h6) 


where 

tn  is  the  time  of  the  nth  sampling  instant, 

x(tn)  is  an  r  X  1  vector  which  is  the  state  of  the  system, 

u(t  )  is  an  r  X  1  vector  which  is  the  effect  of  the  random 
n 

inputs  to  the  system  between  the  n  and  n+1  sample 
points, 

y(t^)  is  a  p  X  1  vector  which  is  the  observed  random  variable, 

$(t  ,,t  )  is  the  r  X  r  transition  matrix, 

v  n+1  n'  ’ 

M(t^)  is  the  r  X  p  output  matrix. 

A  matrix  block  diagram  of  the  sampled  system  is  presented  in  Fig.  k. 

The  element  marked  DELAY  permits  the  state  of  the  system  to  change  only 
at  the  sample  points. 

For  a  discrete  stochastic  process  the  properties  of  the  sampled 
excitation  u(^n)  ^y  be  given  directly,  but  they  can  also  be  derived 
from  the  properties  of  the  random  input  to  the  continuous  system  v(t). 
By  comparing  Eqs.  (4-6)  and  (10)  it  is  seen  that 


u(t  )  =  ( 

n  J 


n+1 


$(tn+1,T)  G(t)  v(t)  dT. 


(V7) 


Thus,  the  sampled  excitation  u('tn)  also  has  zero  mean,  and  the 
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delta -function  property  of  the  covariance  of  the  random  input  v(t)  in 

Eq.  (2)  means  that  the  covariance  matrix  of  the  excitation  is 

t 


t  . 
r  n+l 


E[u(tn)uT(tn)]  =  E 


n+l 


■l(tn+l,T)G(T)v(T)dT 


vT(Cr)GT(o-)ir(tn+1,cr)d( 


n+l 


#  (t  ,,T)G(T)V(t)GT(T)  fT(t  ,T)dT. 


m 


The  sampled  excitation  is  an  independent  random  variable,  so  from 
Eq.  (4) 


E[u(tn)xT(t)]  =  0  tn>t.  (49) 

In  this  investigation,  the  output  matrix  M  and  the  transition 
matrix  f  of  the  model  in  Fig.  4  are  not  known  exactly,  but  they  are 
considered  to  be  matrix-valued  random  parameters  of  the  model  with 
known  probability  distribution  independent  from  one  sample  point  to 
the  next.  The  matrices  M  and  t  will  be  random  parameters  if,  for 
instance,  the  observed  random  variable  has  been  corrupted  by  multipli¬ 
cative  noise  or  if  the  sample  period  is  a  random  variable. 
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More  precisely,  it  is  assumed,  that  a  known  cumulative  distribution 

function  F^[C,n]  is  defined  over  the  components  of  M(tn),  and  a 

different  known  cumulative  distribution  function  Fj[D,n]  is  defined 

over  the  components  of  4  (tn+1,t  ).  When  PR[A]  is  defined  as  the 

probability  that  an  event  A  occurs,  and  B  <  C,  where  B  and  C  are 

matrices,  is  defined  to  imply  that  b. .  <  c3J  for  all  i  and  i,  then 

ij  -  ij 

PR[M(tn)  <  C]  =  PM[C,n] 

PR[#(tn+l,tn)  -  D]  =  (50) 

where  F^[C,n]  and  F^[D,n]  are  the  known  cumulative  distribution 
functions.  These  distribution  functions  are  not  necessarily  stationary, 
but  the  random  parameters  must  be  statistically  independent  from  one 
sampling  instant  to  the  next  so  that 

PR[M(tm)  <  Ai  tj  <  Bj  M(tn)  <  Cj  <  D] 

*  VA’”1  FM[C,n]  F,tD,n).  (51) 

Finally,  it  is  assumed  that  E[M(tn)  M^(tn)],  and  E[4 (tn+^,t  ) 

5*^(t  ^,t  )]  exist  for  all  n. 

The  expected  values  of  the  two  matrices  M  and  4  will  always  exist, 
also,  by  the  previous  assumption,  and  this  expected  value  will  be  denoted 
by  a  superscript  bar. 


E[M(tn)]  =  M(tn) 

-  *  <Wa>-  <52) 

The  difference  between  the  actual  value  of  the  matrix  and  the  expected 
value  will  be  denoted  by  a  superscript  tilde,  so 
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(53) 


M(tn)  =  M(tn)  +  “ftn) 

i  (t  ,t  )  =  s(t  ,,t  )  +}  (t  1  t  ). 

'  n+17  n7  n+17  n7  '  n+1,  n7 

In  this  investigation  it  is  not  necessary  to  know  the  probability 
distribution  of  the  random  parameters.  The  only  information  which  is 
needed  is  the  mean  of  the  parameters,  M(tn)  #(t  ^t^,  and  the 

covariances 


E[M(tn)x(tn)xT(tn)M(tn)] 

and 

E[f(t  .  ,t  )x(t  )xT(t  )  *(t  .  ,t  )]. 

n+17  n  n  v  n7  *  '  n+1  n  J 

E.  EXAMPLES  OF  RANDOM  PARAMETERS 


The  probability  distribution  of  the  random  parameters  M  and  i 
in  the  sampled  model  in  Fig.  k  may  have  any  form,  but  there  are  some 
common  probability  distributions  which  can  be  used  to  approximate  situa¬ 
tions  occurring  in  actual  practice.  In  this  section  the  mean  and  the 
covariance  of  the  random  parameters  will  be  determined  for  five  of  these 
situations;  two  cases  of  multiplicative  noise  (amplifier  noise  and 
exponential  noise),  and  three  cases  of  randomness  in  the  sample  periods 
(periodic  sampling  with  jitter,  periodic  sampling  with  independent 
misses,  and  purely  random  or  "Poisson"  sampling).  For  a  complete  dis¬ 
cussion  of  the  probability  laws  used  to  approximate  these  situations 
see  Parzen  [Ref.  23]- 

In  the  first  two  cases  a  single  state -variable  is  corrupted  by 
multiplicative  noise  and  the  results  can  be  extended  to  the  multivariable 
case.  Let  €  be  an  independent  random  variable  with  normal  distribu¬ 
tion.  The  cumulative  distribution  function  is 


(5b) 


21 


where  f  (v)  is  the  probability-density  function.  The  random  variable 

e  2 
e  has  mean  zero  and  variance  cr  and 

E  [eXe]  =  (55) 

One  source  of  multiplicative  noise  occurs  when  the  state  variable 
passes  through  an  amplifier  with  noisy  gain.  Then 

y  (t  )  =  A(1  +  e)  x  (t  ), 
in  in 

where  A  is  a  constant.  The  mean  and  covariance  are 


E[M(t  )]  =  E  [A(l  +  €)]  =  A 
E[M(tn)x(tn)xT(tn)Mr(tn)]  =  A2  E[e2}  E[x2(tn)] 

=  A2cr^E[x^(tn)  ]  (56) 

Another  source  of  multiplicative  noise  is  the  result  of  taking  the 
exponential  of  a  state  variable  with  added  white  noise.  In  a  communica¬ 
tion  system  in  which  the  range  of  the  state  variable  is  several  orders 
of  magnitude  (such  as  the  reflected  pulse  in  a  radar  system  where  the 
energy  is  inversely  proportional  to  the  fourth  power  of  the  distance), 
it  may  be  desirable  to  transmit  the  logarithm  of  the  state  variable  and 
then  convert  this  at  the  receiver  to  the  estimate  of  the  original  state 
variable.  In  that  case, 


log  y  (t  )  =  log  x  (t  )  +  e 
in  in 


yL(tn)  =  e£  x^tj 


(57) 


and  the  mean  and  covariance  are 


2? 


(58) 


E[M(tn)]  =  E[e£] 
E[M(tn)x(tn)xT(tn)Mr(tn)l 


=  //* 

=  {jS[e2€]  -  (E[e€l^|  E[x2(tn)] 
=  e*  (e*  -  1)  E[x2(tn)] 


In  the  remaining  three  cases,  the  transition  matrix  (t  . ,t  )  is 
the  random  parameter  because  the  sample  period  is  an  independent 

random  variable.  It  will  be  assumed  that  the  dynamics  matrix  F  is 
constant  or  can  be  satisfactorily  approximated  by  a  constant  over  the 
sample  period  T  ,  but  the  constant  and  the  probability  distribution  of 
T  can  change  from  one  sample  point  to  the  next,  so  the  sample  period 
is  not  necessarily  stationary.  The  transition  matrix  (t  +p,tn)  is 
an  r  X  r  matrix-valued  random  variable  given  by 


FT 


(t  , ,t  )  =  e 
v  n+1'  n' 


(59) 


The  exponential  of  the  matrix  F  T^  is  composed  of  terms  of  the  form 


X.T  X.T  0  X.T  X.T 

i  n  m  i  n  m2  in  _v  in 

e  ,Te  ,Te  ,...,Te 
’  n  n  n 

where  is  a  characteristic  root  of  F  of  multiplicity  v.  The 

powers  of  T  occur  when  X,  is  a  multiple  characteristic  root,  as  in 
n  l 

the  2X2  matrix  in  Eq.  (33)’ 

In  order  to  determine  the  expected  value  and  the  covariance  of  the 
transition  matrix,  i+  is  necessary  only  to  know 


XT 

E[e  n] 


for  all  values  of  X  both  real  and  complex.  If  that  expected  value  is 
known,  the  other  expected  values  can  be  calculated  by  the  relationship 


X.T  d  XT 

E[T  V  e  1  n]  =  —  E[e  n] 
dX 


X  =  X. 
1 


(6o) 
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The  three  situations  in  which  the  sample  period  can  be  considered 
as  a  random  variable,  and  their  associated  probability  distributions, 
are  as  follows:  periodic  sampling  with  jitter  (normal  distribution), 
periodic  sampling  with  independent  misses  (geometric),  and  purely 
random  or  Poisson  sampling  (exponential).  The  pertinent  information 
concerning  these  three  probability  laws  is  presented  in  Table  1,  and 
their  characteristics  are  discussed  in  the  remainder  of  this  section. 


TABLE  1.  THREE  PROBABILITY  LAWS  FOR  RANDOM  SAMPLE  PERIODS. 


Probabi  li  ty 

Law 

Normal 

Geometri c 

Expon  en  ti  al 

Sampl e 

Pattern 

Periodic  with 
Jitter 

Peri odi  c  wi  th 
Indep endent 
Misses 

Purely  Random 

or 

Poi sson 

Parameters 
of  Probability 

Di stri bu  tion 

0  <  30-  <  T  <  “ 

0  <  T  <  00 

0  <  p  =  1-q  <  1 

0  <  yC  <  CO 

Probabi  li  ty 

Den si  ty 
f  (T) 

exp  [-(T  -  T  )  2/20-2] 
/2n  a 

pqk*  *  for  T  =  kT 
k  =  0,  1,  2..  .  . 

0  Otherwise 

Me'^7-  for  T  >_0 

0  Otherwise 

Mean 

E[t] 

T 

T/p 

l/M 

Vari an  ce 

E[t2]  -  £2[t] 

0-2 

qT2/p2 

1/m2 

E[eXT] 

eT\e*A2 

peXT/(l  -qe^T) 

A/Cm-M 

Jitter  occurs  when  the  sample  period  is  neatly  constant,  but  varies 
slightly  from  period  to  period.  In  certain  anti -jamming  applications 
the  sample  period  is  varied  randomly  in  this  manner,  or  the  variation 
may  be  unintentional  due  to  imperfections  in  the  equipment.  An  approxi¬ 
mation  to  the  effect  of  jitter  is  given  by  the  normal  distribution  of 
sample  periods. 


2  2 

F  [t]  =  T  +  -  j  e"V  /2o  dv  0  <  3o  <  T  <  oo  (6l) 

n  f2n  a  ^ 
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or,  alternatively, 


Tn  =  T  +  €,  (62) 

where  T  is  the  mean  of  the  sample  period,  and  e  is  the  normal  random 
variable  with  mean  zero  and  variance  defined  in  Eg.  (5*0*  In  the  normal 
distribution  in  Eg.  (6l),  there  is  a  finite  probability  that  the  sample 
period  will  be  negative;  therefore,  one  might  argue  that  a  better 
approximation  would  be  to  truncate  the  distribution  so  that  it  is  zero 
outside  the  region 


0  <  T  <  2T. 
n 

This  is  a  valid  argument,  but  in  Eg.  (6l)  the  standard  deviation  is 
constrained  to  be  less  than  one-third  of  the  mean  of  the  sample  period 
T  and,  under  this  constraint,  the  area  of  the  normal  distribution  that 

_  O 

must  be  truncated  is  less  than  2.6  X  10  therefore,  the  original 
approximation  should  be  valid. 

Periodic  sampling  with  independent  misses  results  when  the  intended 
sample  period  is  constant  but  at  each  sampling  instant  there  is  a  fixed 
probability  g  that  the  random  variable  will  not  be  observed.  This 
situation  may  arise  when  the  receiver  rejects  the  signal  unless  the 
signal-to-noise  ratio  is  above  some  threshold  value.  The  probability 
distribution  of  the  actual  sample  period  is  a  geometric  distribution 
given  by 


fT  (kT)  =  pg(K_x)  k=l,2,...  0<p=l-g<l 

n 

=  0  otnerwise  (63) 

When  the  number  of  sample  points  in  a  given  interval  of  time  has  a 
Poisson  distribution,  the  samples  are  being  received  at  purely  random 
points  in  time.  The  samples  may  be  transmitted  in  this  way  intentionally 
to  avoid  jamming,  or  because  of  the  random  character  of  the  guantity 
being  transmitted.  For  this  case  the  sample  period  has  an  exponential 
distribution 
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(6k) 


Ft  (t)  =  n  e"^T  T  >  0  0  <  n  <  « 

n  =0  otherwise 

F.  THE  OFTIMJM  ESTIMATE 

The  purpose  of  the  design  procedure  presented  in  this  paper  is  to 
estimate  the  state  vector  of  a  stochastic  process  from  information 
obtained  from  the  a  priori  knowledge  of  the  process  and  the  observed 
random  variables.  The  difference  between  the  actual  value  of  the 
state  vector  and  the  estimate  is  the  error,  which  is  written 

x(tQ  +  a/tn)  =  x(tn  +  a)  -  x(tn  +  a/tn),  (65) 

where 

x(t  +  a/t  )  is  the  (r  X  l)  vector  that  is  the  error  in  the 
estimate  of  the  value  of  the  state  vector  at  time 
t^  +  a,  given  the  observed  random  variables 
y(tn),y(tn_1), . . .,  y(tQ); 

x(t  +  a)  is  the  actual  value  of  the  state  vector, 

~  n 

x(t  +  a/t  )  is  the  estimate  of  the  value  of  the  state  vector 
n  n 

at  time  t  +  a,  given  the  observed  random  variables 
y(tn),  y(tn_1),...,  y(tQ). 

The  best  estimate  is  the  one  which  minimizes  some  function  of  the  error. 
For  a  positive  the  operation  of  estimation  is  called  prediction;  and 
for  a  negative,  it  is  called  interpolation  or  smoothing. 

In  this  paper  the  estimate  will  be  confined  to  linear  combinations 
of  the  observed  random  variable,  and  the  optimum  estimate  will  be  the 
linear  estimate  that  minimizes  the  trace  of  the  covariance  matrix  of 
the  error.  The  covariance  of  the  error  is  an  r  X  r  matrix  defined  as 

p(tn  +  a/\)  =  E[x(tn  +  a/tn)xT(tn  +  a/tn)].  (66) 

The  trace  is  the  sum  of  the  diagonal  terms  of  the  matrix  The  trace 
of  the  covariance  is  the  expected  value  of  the  sum  of  the  squared  error 
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in  the  estimation  of  the  individual  state  variables,  so  the  optimum 
estimate  is  a  natural  extension  of  the  "linear-least-squares"  estimate 
that  is  discussed  extensively  in  the  literature. 

The  question  naturally  arises  as  to  which  is  "better — a  linear  esti¬ 
mate  or  a  nonlinear  estimate.  Dooh  [Ref.  24]  proves  rigorously  that  the 
estimate  of  a  quantity  that  minimizes  the  expected  value  of  the  squared 
error  is  the  conditional  expectation  of  the  quantity  with  respect  to  the 
given  information.  For  a  linear  dynamic  system  excited  "by  random 
variables  with  normal  distribution,  the  conditional  expectation  of  the 
state  vector  is  indeed  a  linear  combination  of  the  observed  random 
variables  and  the  linear-least-squares  estimate  will  be  better  than  any 
nonlinear  estimate  in  minimizing  the  expected  value  of  the  squared  error. 

Sherman  [Ref.  3]  has  proved  that  the  conditional  expectation  can 
minimize  the  expected  value  of  a  still  larger  class  of  loss  functions. 

In  particular,  when  e  is  the  random  variable  representing  the  error, 
and  the  loss  function  L(e)  is  a  positive  function,  symmetric  and  non¬ 
decreasing  about  zero,  so  that 

0  =  L(e)  =  L(-e) 

L(e1)  =  L(e2)  when  0  =  ex  =  e2,  (67) 

and  if  the  conditional  distribution  of  the  quantity  being  estimated  is 
symmetric  and  convex  about  its  conditional  expectation,  then  the  condi¬ 
tional  expectation  minimizes  the  expected  value  of  the  class  of  loss 
functions  included  in  Eq»  (67)  This  class  of  loss  functions  includes: 

Lx(e)  =  e2 

L2(e)  =  e^ 

L3(e)  =  |e | 

L^(e)  =1  -  b  <  e  <  b 

=  0  otherwise. 
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(68) 


In  this  investigation,  Sherman's  results  mean  that  if  the  stochastic 
process  is  the  output  of  a  linear  dynamic  system  with  a  gaussian  input, 
and  if  the  random  parameters  have  a  symmetric  convex  probability  distri¬ 
bution,  then  the  optimum  linear  estimate  will  be  better  than  any  other 
estimate,  linear  or  nonlinear,  in  minimizing  the  expected  value  of  the 
class  of  loss  functions  in  Eq.  (67). 

The  optimum  linear  estimate  has  the  additional  advantage  that  only 
the  first-  and  second-order  statistical  averages  are  required.  In 
practical  situations  it  is  often  very  difficult  to  measure  more  than 
this,  and  any  nonlinear  estimate  would  require  more  information  than 
just  the  first-  and  second-order  statistical  averages  to  improve  over 
the  optimum  linear  estimate.  For  these  reasons,  in  this  investigation 
the  optimum  estimate  will  be  the  linear  estimate  that  minimizes  the 
trace  of  the  covariance  matrix  of  the  error. 
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III.  OPTIMUM  FILTERING  AND  PREDICTION  OF  RANDOM- PARAMETER  PROCESSES 


A.  INTRODUCTION 

In  Chapter  II  the  sampled  stochastic  process  was  represented  as  a 
matrix-difference  equation  with  random  parameters.  The  optimum  estimate 
was  defined  as  the  linear  estimate  that  minimized  the  trace  of  the  co- 
variance  matrix  of  the  error.  In  this  chapter  it  will  be  proved  that 
the  optimum  estimate  of  the  current  value  of  the  state  vector  at  the 
current  sample  point  is  a  linear  combination  of  the  optimum  estimate  at 
the  preceding  sample  point  and  the  random  variable  observed  at  the 
current  sample  point.  The  optimum  estimate  is  implemented  by  a  linear 
filter  with  the  only  undetermined  coefficient  being  a  matrix- valued 
gain.  This  matrix- valued  gain  as  well  as  the  covariance  of  the  error 
in  the  optimum  estimate  is  determined  by  a  nonlinear  difference  equation 
that  can  be  solved  iteratively  for  each  sample  point.  The  optimum 
current  estimate  becomes  the  optimum  prediction  of  a  future  value  of 
the  state  vector  when  it  is  matrix-multiplied  by  the  expected  value  of 
the  transition  matrix. 

These  results  are  obtained  from  the  solution  to  the  matrix-valued 
sampled  Wiener-Hopf  equation,  derived  in  Sec.  III-C,  but  first  the 
Markov  property  of  the  optimum  linear  estimate  will  be  discussed. 

B.  MARKOV  PROPERTY  OF  THE  OPTIMUM  ESTIMATE 

For  a  Markov  process,  the  probability  functions  relating  to  the 
future  depend  on  the  present  state,  but  not  on  the  manner  in  which  the 
present  state  has  emerged  from  the  past  [Ref.  25],  The  model  of  the 
stochastic  process  developed  in  Chapter  III  has  the  Markov  property, 
and  in  Sec.  III-D  it  will  be  proved  that  the  optimum  linear  estimate 
has  the  Markov  property  because  the  optimum  estimate  of  the  current 
value  of  the  state  vector  at  the  current  sample  point  is  a  linear  combi¬ 
nation  of  the  optimum  estimate  ax  the  preceding  sample  point  and  the 
observed  random  variable  at  the  current  sample  point  In  other  words, 
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/ 


^VV  =  *  *(tn^tn-l)$(tnVtn)  +  K(tn)y(tn)  (69) 

_  * 

where  <t>  (tn,tn_^)  and  K(tn)  are  r  x  r  and  r  x  p  linear 

weighting  matrices.  (Note  that  0  (t  ,t  is  an  arbitrary 

weighting  coefficient  different  from  the  expected  value  of  the  transi¬ 
tion  matrix  ‘t>  (t  ,t  ,  ),  which  does  not  have  the  asterisk.)  This  see- 
n  n-1  7 

tion  will  explore  the  reasoning  which  might  lead  one  to  the  hypothesis 
of  Eq.  (69). 

First,  assume  that  at  time  t  .  the  exact  value  of  the  state 

n-1 

vector  is  known;  therefore, 

x(t  /t  )  =  x(t  ).  (70) 

n-1  n-1  n-1  ' 

If  the  transition  matrix  in  the  sample  period  is  equal  to  its  expected 
value  $  (t^t^  -^)  and  if  there  are  no  random  inputs  to  the  system 
during  the  sample  period,  then  the  next  state  will  be 

x(tn)  =  *  ( W^xtt^),  (71) 

and  the  best  estimate 

x(t  )  =  4>  (t  ,t  )x(t  )  =  x(t  )  (72) 

n  '  n7  n-1  n-17  n7  '  ' 

will  be  the  exact  value  of  the  state  vector.  If  the  output  matrix  is 

also  equal  to  its  expected  value  M(t  ),  the  observed  random  variable 

n 

will  be 


y(t.  )  =  M(t  )x(t  )  =  M(t  )  0  (t  ,t  )x(t  ,).  (73) 

and  no  new  information  about  the  process  will  be  gained  by  looking  at 
it,  because  its  exact  value  is  already  known.  This  is  not  a  very 
interesting  problem,  and,  naturally,  in  actuality  things  are  not  so 
simple.  There  will  be  errors  in  the  estimate  of  the  state  at  the  pre¬ 
various  sample  instant,  x(t  , /t  .  );  there  are  random  excitations  to 

n-1'  n-1 

the  system  u(t^  and  variations  in  the  transition  and  output  matrices 

$  (t  ,t  n)  and  M(t  ). 

'  n  n-1  '  n 
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The  only  way  one  discovers  these  random  variations  is  by  comparing 
the  observed  random  variable  y(t  )  with  its  expected  value  which  is 
given  by  Eq.  (73).  All  the  new  information  about  the  process  comes  from 
the  difference  between  y(t  )  and  its  expected  value 

A 

y(t  )  -  M(t  )  $  (t  ,t  . )x(t 

n'  n'  v  n  n-1'  '  n-17 

Because  the  estimate  is  restricted  to  linear  operations,  it  is  reasonable 
to  weight  ^he  new  information  with  some  linear  matrix-valued  coefficient 
K(tn)  and  add  it  to  the  previous  estimate  in  Eq.  (72),  so  that 

x(t  /t  )  =  <J>  (t  /t  -i  )x(t  /t  ) 

'  n'  n  '  n'  n-1'  v  n-1'  n-1 

+  K(t  )[y(t  )  -  M(t  )  4>  (t  /t  )x(t  /t  )}.  (7*0 

v  n'L,/'  n'  '  n'  n'  n-1'  '  n-1'  n-1  J  ' 


Combining  terms  in  Eq.  (7*0  gives 

x(t  /t  )  =  [I  -  K(t  )M(t  )]  $  (t  ,t  -i  )x ( t  ,/t  )  +  K(t  )y.(t  ) 

'  n'  n'  1  '  n'  v  n'J  v  n'  n-1'  v  n-1'  n-1  n  v  n' 

where  I  is  the  identity  matrix  and 

[I  -  K(tn)M(tn)]  9  (tn,tn_1)  =  0  \ta,tn_1),  (75) 

the  arbitrary  r  x  r  weighting  matrix  in  Eq.  (69). 

This  reasoning  started  with  the  assumption,  in  Eq.  (70),  that  at 

time  t  ,  the  exact  value  of  the  state  vector  was  known,  but  all  the 
n-1  ' 

ideas  still  apply  if  only  the  expected  value  of  the  state  vector  is 

known.  The  key  idea  in  this  discussion  can  be  summarized  as  follows: 

When  the  stochastic  process  can  be  represented  as  an  r^-order  linear 

Markov  process,  the  optimum  linear  estimate  of  the  state  vector  of  the 

th 

process  can  also  be  represented  as  an  r  -order  linear  Markov  process. 

This  approach  to  the  linear  estimate  of  a  stochastic  process  was 
first  developed  fully  by  Kalman  [Ref.  1]  for  a  process  with  deterministic 
parameters.  He  used  conditional  expectations  and  the  projection  theorem, 
and  considered  in  detail  the  linear  estimate 
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(76) 


$(tn+l/tn)  =  *  +  AX)y(V’ 


where 


*. 

$  (t 


,t  )  =  $ (t  ,  ,t  ) 
n+1  n'  v  n+1'  n' 


*. 

A  (t 


)M(t  ). 
n  v  n' 


(77) 


He  proved  that  any  other  estimate  x(tn  +  a/t^)  for  a  positive  could 
be  derived  from  the  one  in  Eg.  (76)  by  the  relationship 

x(t  +a/t  )  =  <t(t  +a,t  )[®(t  ,  ,t  )]  'Lx(t  .,/t  ),  (78) 

'  n  '  n'  '  n  '  n/l  v  n+x’  n/J  v  n+1'  nJ>  v  ‘ 

where  A  represents  the  inverse  of  A. 

The  quantity  that  in  this  investigation,  is  analogous  to  the  quan¬ 
tity  in  brackets  in  Eq.  (78)  is  (tn+^,tn),  which  may  be  singular 
and  does  not  always  have  a  unique  inverse.  Therefore,  for  the  stochastic 
process  with  random  parameters  this  approach  must  be  modified  to  the 
linear  estimate  in  Eq.  (69). 

C .  THE  SAMPLED  WIENER-HOPF  EQUATION 

Any  linear  estimate  of  the  state  vector  will  be  a  linear  combination 
of  the  observed  random  variables  so  that 

n 

*(tn+a/tn)  =  A(tn+a'tv)y(tv)>  (79) 

v=0 

is  an  r  x  1  vector  that  is  the  estimate  of  the 
state  vector  x(tn+a),  given  the  observed  random 
variables  y(tn),y(tn_1), . . . ,y(tQ  ),  and 
is  an  r  x  p  matrix  which  is  the  vth  set  of 
weighting  coefficients  of  the  estimate. 

The  sampled  Wiener-Hopf  equation  is  a  matrix- valued  linear  equation 
that  is  satisfied  by  the  weighting  coefficients  A(tn+a,t^)  when  the 
linear  estimate  is  optimum. 


where  x(t  +a/t  ) 
'  n  '  n' 


A(t  +a,t  ) 
'  n  v 
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The  error  in  the  linear  estimate  in  Eq.  (79),  as  defined  in  Eq.  (65), 


is 


n 

X(tn+a/tn^  =  x(tn+a)  *  x(tn+a/tn)  =  x(tn+a)  -  ^  A(tn+a,tv)y(tv) 

v=0 

(8o) 

Substituting  the  error  in  Eq.  (80)  into  the  covariance  of  the  error 
defined  in  Eq.  (66)  gives 

p(tn+a/tn)  =  E[x(tn+a/tn)xT(tn+a/tn)] 

=  E[x(tn+a)xT(tn+a)] 

n 

-T  A(tn+a,tv)E[y(tv)xT(tn+0()] 

v=0 

n 

E[x(tn+a)yT(tv)]AT(tn-fO!,tv) 

v=0 
n  n 

+  y  y  A(tn+a^tv)  E[y(tv)yT(t^)]  ^(t^,^)  (8l) 

v=0  n=0 


The  trace  of  the  covariance  matrix  in  Eq.  (8l)  is  the  sum  of  the 
diagonal  terms,  so  that 


r 

TR[P(tn+<Vtn)]  =y  Pii(tn+a/tn)  =  E 

i=l 

which  is  the  expected  value  of  the  sum  of  the  squared  error  in  the 

estimation  of  the  individual  state  variables.  From  Eq.  (8l)  it  is  seen 

that  all  the  components  of  the  covariance  matrix  are  quadratic  functions 

of  the  elements  of  the  weighting  coefficients  a„  (t  +a, t  ),  so  the 

jk'  n  v 

trace  in  Eq.  (82)  is  minimized  when  its  derivative  with  respect  to  each 

of  the  elements  a.,  (t  +a,t  )  is  zero, 
jk'  n  v 


X 

y 


i=l 


X^(t  *Jt/t  ) 


(82) 
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(83) 


d  TR[P(tnK*/tn)] 
d  a.k(tn40!,tv) 


j  -  1,2, ...  ,r 
k  =  1,2, ... ,p 
v  =  0, 1, . . . ,n 


Because  the  trace  is  a  quadratic  function  of  the  elements,  there  will 
be  at  least  one  minimum,  and  the  resulting  equation  for  the  optimum  ele¬ 
ments  will  be  linear.  After  performing  the  operation  indicated  in  Eq. 
(83)  and  combining  the  linear  equations,  the  result  is 

n 

E[x(tn+a)yT(tp)]  -  ^A(tn+a,tv)E[y(tv)yT(tp)]  =0  p  =  0,1,  ...,n, 
V=0  (84) 

which  is  the  sampled  Wiener-Hopf  equation.  This  result  yields  a 
minimum  because 


d  TR  [  P  ( t  +0t/t  )  ] 

- - V-E[y2(t  )]>0 

d  ajk(tn40,tv) 


(85) 


the  second  derivative  with  respect  to  the  elements  is  positive.  The 
optimum  weighting  in  Eq..  (81+)  not  only  minimizes  the  trace,  but  it  also 
minimizes  all  the  components  of  the  covariance  matrix,  as  can  be  seen  by 
carrying  through  the  operations  in  Eq.  (83)  for  the  off-diagonal  compo¬ 
nents  of  the  covariance  matrix.  Substituting  the  estimate  in  Eq.  (79) 
for  the  equivalent  expression  in  Eq.  (84)  gives 


E[x(tn+a)yT(tp)]  -  E[x(tn+a/tn)yT(tp)]  =  0 


p  =  0, 1, . . . ,n  (86) 


as  a  compact  way  of  writing  the  sampled  Wiener-Hopf  equation.  When  the 
weighting  coefficients  are  optimum  and  satisfy  Eq  (86),  some  of  the 
terms  in  the  expression  for  the  covariance  matrix  in  Eq.  (8l)  are 
identically  zero,  and  the  covariance  of  the  error  can  be  written, 


P(tn+a/tn)  =  E[x(tn+a)xT(tn+a)]  -  E[x(tn+a/tn)xT(tn-+o:)] 


where  x(tn+a/tn)  i-s  defined  in  Eq.  (79)- 


(87) 
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D.  DERIVATION  OF  THE  OPTIMUM  ESTIMATE 


The  sampled  Wiener-Hopf  Eq.  (86)  must  be  satisfied  by  the  weighting 

coefficients  of  the  optimum  linear  estimate  for  any  stochastic  process 

with  any  statistics.  In  this  section  it  will  be  proved  that,  when  the 

process  can  be  represented  by  the  model  developed  in  Chapter  II,  the 

optimum  estimate  has  the  Markov  property  discussed  in  Sec.  III-B.  The 

assumption  is  first  made  that  the  optimum  estimate  of  the  current  value 

of  the  state  vector  x(t  /t  )  can  be  written  in  the  form 

n  n 

x(t  /t  )  =  [I  -  K(t  )M(t  )]  j>  (t  ,t  .  )x(t  /t  )  +  K(t  )y(t  ), 
n'  n'  1  '  n  x  n  J  x  n  n-1  v  n-1'  n-1'  n  n" 

(88) 

where  x(t^  ^/t  ^)  is  the  previous  optimum  estimate  and  K(tn)  is  an 

r  x  p  weighting  matrix.  Next  it  is  shown  that,  if  the  weighting  coef¬ 
ficients  of  the  previous  estimate  x(tn  ^An  ^)  satisfy  the  sampled 
Wiener-Hopf  Eq.  (86),  the  weighting  of  the  current  estimate  in  Eq  (88) 
will  satisfy  Eq.  (86)  for  a  particular  r  x  p  weighting  matrix  K(t^). 
Finally,  the  weighting  matrix  K(tn)  is  derived  in  terms  of  a  priori 
known  quantities  and  the  covariance  of  the  error  of  the  previous  optimum 
prediction  P(tn/tn  ^),  which  can  be  calculated  iteratively  at  each 
sample  point. 

The  random  inputs  to  the  system  have  zero  mean;  therefore,  from 
Eq.  (ll)  for  a  positive, 

E[x(tn+a)yT(tp)]  =  E[4>(tn+a,tn)x(tn)yT(tp)]=  $(tn+a,tn)E[ x(tn)yT(tp) ] . 

a  >  0 

P  <  n  (89) 

Substituting  Eq.  (89)  into  Eq.  (86)  for  a  positive  gives 

o(t  +a,t  )Efx(t  )y'1(t  )]  -  E[x(t  +a/t  )yT(t  )]  =  0  a  >  0 
'  n  n  L  n  p  J  1  '  n  '  n  p  1  — 

p  =  0,1, . . .n 

(90) 
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as  the  condition  for  optimum  prediction  It  is  obvious  from  Eq.  (90) 

that,  if  x(t  /t  )  is  the  optimum  current  estimate  and  hence  satisfies 
'  n'  n 

that  equation  for  a  equal  to  zero, 

x(t  +a/t  )  =  $  (t  +a,t  )x(t  /t  )  a  >  0  (91) 

n  n  n  n  v  n'  n  -  x  7 

must  satisfy  Eq.  (90)  for  a  positive.  Thus,  the  optimum  prediction 
is  obtained  from  the  optimum  current  estimate  by  matrix  multiplying  the 
current  estimate  by  the  expected  value  of  the  transition  matrix.  Now, 
the  optimum  current  estimate  will  be  determined. 

Assume  the  previous  optimum  estimate  is 


n-1 

x(t  /t  )  =  /  A(t  ,t  )y(t  ),  (92) 

'  n-1'  n-17  i_j  '  n-17  v7  v  7 

0 

and  therefore  it  satisfies  Eq.  (93  )> 

Etx(tn.i)yT;tp)]  -  E^(tn-l//tn-l^yT(tp^  “  0  p  =  0,1,..., n-1,  . 

(93) 

From  Eqs.  (91 )  and  (92),  the  optimum  prediction  at  the  previous  sample 
instant  of  the  current  value  of  the  state  vector  is 

x(t  /t  )  =  <J>  (t  ,,t  )x(t  ,/t  .  ),  (94) 

v  n'  n-1  v  n+1  n'  v  n-1'  n-i  v 

and 'from  Eq.  (87)  the  covariance  of  the  error  of  this  prediction  is 

P(t  /t  )  =  E[x(t  )xT(t  )]  -  d>(t  ,t  )E[x(t  /t  )xT(t  )]. 
n  n-1  n'  n,J  '  n  n-1  n-1'  n-1  v  n'J 


The  trial  estimate  in  Eq  (88)  can  be  written 


x(t  /t  )  =  (t  ,t  .  )x(t  /t  ) 
n  n  n  n-1'  v  n-1'  n-1 


+  K(tn)[y(tn)  -  M(tn)(t  ^n’Vl^^n-A-l”1 


(95) 


(96) 
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A 

For  the  estimate  x('tnAn)  Eq.  (96)  to  be  optimum,  it  must  satisfy 
Eq.  (97), 

E[x(tn)yT(tp)]  -  E[x(tn/tn)yT(tp)]  =  0  .  p  =  0, 1, . . . ,n.  (97) 

Substituting  the  trial  estimate  in  Eq.  (96)  into  Eq.  (97)  gives 

E[x(tn)yT(tp)]  -  $  (Vtn-l)E^(tn-l/tn-l)yT(tp)] 

-  K(tn)^[y(tn)yT(tp)]  -  M(tnJ  ^Vtn_1)E[x(tn_1/tn_1)yT(tp)]J  =  0 

p  =  0,1, .,n  (98) 

as  the  condition  for  optimum  estimation. 

Equation  (98)  will  be  examined  first  for  p  less  than  n.  In  that 
case, 

E[x(t  )yT(t  )]  =  $(t  ,t  )E[x(t  )yT(t  )] 

1  p<J  n  n-1  n-1  p  J 

E[y(tn)yT(tp)]  =  M(tn)  ^(tn,tn_i)E[x(tn_i)yT(tp)]  P  <  n 

(99) 

Substituting  Eq.  (99)  into  Eq.  (98)  for  p  less  than  n  yields, 

[I  -  K<tn)M(tn)]  ^('tn^tn_1)^:[x(tn.1)yT(tp)]  -  E[x(tn.i/tn-l^T(tp)l|=  0 

p  =  0,1, ...,n-l  (lOO) 

Comparison  of  Eqs.  (100)  and  (93)  shows  that  the  quantity  in  braces  in 
Eq.  (lOO)  is  the  same  as  Eq,  (93)>  which  is  identically  zero  for  p 
less  than  n.  Thus,  the  trial  estimate  in  Eq.,  (95)  satisfies  the  condi¬ 
tions  for  optimum  estimation  in  Eq,  (98)  for  p  less  than  n.  The 
remaining  condition  in  Eq.  (98)  for  p  equal  to  n  will  be  used  to 
determine  the  r  x  p  weighting  matrix  K(t^);  this  condition  is 
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EMtn)yT(t„)]  -  »(tn.tn.1)s[5(t„.1/tn..1)yT(tn)] 

-  K(tn)^[y(tn)yT(tn)]  . 

(101) 

From  the  model  of  the  process  in  Eq.  (46)  and  the  properties  of  the 
random  parameters  in  Eq.  (53)  *  the  observed  random  variable  y(t  )  can 
be  written 


y(tn)  =  M(tn)x(tn)  =  M(tn)x(tn)  +  ft(tn)x(tn)y  (102) 

where  M('tn)  is  an  independent  random  variable  with  zero  mean,  so  that 

E[x(t  )yT(t  )]  =  E[x(t  )xT(t  )]  M^t  ) 

1  n  n  J  1  N  n  v  n  J  x  n 

E[x(t  /t  )yT(t  )]  =  E[x(t  /t  )xT(t  )]  M^t  ) 

1  v  n-1  n-1  v  n  J  1  v  n-1'  n-1  v  n  J  '  n 

E[y(tn)yT(tn)]  =  M(tn)E[x(tn)xT(tn)]MT(tn) 

+  E[M(tn)x(tn)xT(tn)MT(tn)] 


=  M(tri)E[x(tn)xT(tn)]wr(tn)  +  R(tn). 


(103) 


In  the  third  equation  of  Eq  (103),  the  substitution  has  been  made  that 


R(tn)  =  E[M(tn)x(tn)*T(tn)MT(tn)]. 


(104) 


Substituting  Eq.  (103)  into  the  remaining  condition  for  optimum  estima¬ 
tion  in  Eq.  (lOl)  and  combining  terms  yields 


[I 


-K(tn)M(tn)]|E[x(tn)xT(tn)]-4’(t:i,tn_1)Erx(tn_1/tn_1)xT(tn)]jMr(tn) 


-  K(t  )  R(t  )  =  0 
n'  n 


(105) 


The  quantity  in  braces  in  Eq.  (105)  is  equal  to  Eq.  (95)  which  is  the 
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covariance  of  the  error  in  the  optimum,  prediction  at  the  previous  sample 
instant.  Replacing  that  quantity  by  P(t  /tn  in  Eq.  (95)  yields 

p(t  /t  )if(t  )  -  K(t  )|k(t  )p(t  /t  jifct  )  +  R(t  n  =  o, 
n'  n-1  '  n  n'J  v  n  n  n-1  n  n'J 

(106) 

where  R(tn)  is  defined  in  Eq.  (l04)  When  the  r  x  p  weighting  matrix 
K(t^)  satisfies  Eq,  (106),  the  trial  estimate  in  Eq,  (96)  is  indeed  the 
optimum  estimate  because  it  satisfies  the  sampled  Wiener-Hopf  Eq.  (97). 
The  only  quantity  in  Eq.  (106)  that  is  not  known  a  priori  is  the  covar¬ 
iance  matrix  of  the  error  p(t  /t  0°  The  covariance  P(t  .,/t  )  can 

'  n'  n-1  n+1'  n' 

be  determined  iteratively  at  each  sample  point  from  the  covariance 
P(t  /t  ^)  at  the  previous  sample  point. 

From  Eq.  (9l)  the  optimum  prediction  x(t  ,/t^) 


x(t  ,/t  )  =  ^(t  « ,  t  )x(t  /t  ), 

'  n+1  n'  v  n+1'  n'  v  n'  n" 


(107) 


From  Eqs.  (87)  and  (89)  the  covariance  of  the  error  in  this  prediction 
is 


P(tn+l/tn)  =  Etx^t-xi 


n+1'  '  n+1 


n+1'  n' 


n+1' 


=  E[x(t  n  )xJ'(t  n)l  -  <j>(t  n,t  )ETx(t  /t  )xT(t  )]i^'(t  ..  t  ). 

1  v  n+1'  '  n+1' J  '  n+1  n'  v  n'  n  v  n'J  v  n+1,  n' 

(108) 

The  form  of  the  optimum  estimate  x(tn/t^)  is  given  by  Sq.  (96)  and, 
using  Eq.  (9*0,  this  can  be  written 

x(t  /t  )  =  0(t  ,t  )x(t  /t  )+K(t  )[y(t  )-M(t  )4>(t  ,t  )x(t  /t  )] 
'  n  n  '  n  n-1'  '  n-1  n-1'  n'L  '  n'  '  n  '  n  n-1'  n-1  n-1  1 

=  K(tn)y(tn)  +  [I  -  K(tn)M(tn)]x(tn/tn_1).  (109) 


Substituting  Eq.  (109)  into  the  covariance  in  Eq.  (108)  gives 
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IfcMgg; 


(no) 


From  the  model  of  the  process  in  Eq.  (46)  and  the  properties  of  the 
random  parameters  in  Eq,  (53),  the  state  vector  x(tr+1)  can  he  written 

x(t  )  =  <£  (t  ,t  )x(t  )  +  u(t  ) 

'  n+1'  '  n+1'  n'  '  n'  '  n 

=  $  (t  ,t  )x(t  )  +  $(t  ,t  )x(t  )  +  u(t  ).  (ill) 

n+1  n  n  n+1  n'  '  n'  '  n'  v  ' 

Therefore,  the  covariance  of  x(t  , )  is 
*  n+1 

E[x(tn+l)x^(tn+l)  ]  =  ^Vl'^^n^X^  ^n+l'V  + 


(112) 

where  the  substitution  has  been  made  that 

Q(tn)  =  E[$(tn+1,tn)x(tn)xT(tn)  t(tn+1,tn)]  +  E[u(tn)uT(tn)], 

(113) 

and  the  cross  terms  in  Eq.  (112)  are  zero  because  ^(t^-^t^ )  and  u(t^  ) 
have  zero  mean  and  are  independent  of  x(t^). 

Substituting  Eq.  (112)  into  the  covariance  in  Eq.  (110)  gives 

‘■(WV  ■  ‘  x 

^E[x(t[i)xT(tti)’-ErJ(tn/tn_1)xT(tn)]'j»T(tijtl,tn)  *  5(tn)  (nit) 

The  quantity  in  braces  in  Eq,  (ll4)  is  the  covariance  of  the  error 
p(tn/tn  ■]_)  in  Eq-  (95)-  Replacing  that  quantity  in  Eq.  (ill)  by 
p(t  /t^  ^)  yields  as  the  iterative  equation  for  the  covariance  of  the 


error 
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P<tn.l’t„>  ■  ♦<tn+l'tn)(I  -  K<t„)M(t„>}  P<Vt„-l)  ^‘n+l’V  +  “<*■»>' 

(115) 

•where  Q(tn)  which  is  defined  in  Eq.  (H3)  and  the  two  expected  values 
M(tn)  and  $(tn+^,tn)  are  known  a  priori.  The  covariance  of  the  error 
for  the  current  optimum  prediction  is  given  ty  Eq.  (115)  when  P(t  /t  ) 
is  the  covariance  of  the  previous  optimum  prediction  and  the  r  x  p 
weighting  matrix  K(tn)  satisfies  Eq.  (106)  [which  is  repeated  as 
Eq.  (il6)]. 

P(t  /t  jM^t  )  -  K(t  )  <  M(t  )P(t  /t  Jffct  )  +  R(t  )l  =  0. 

'  n'  n-1'  '  n'  '  n'  1  '  n'  '  n'  n-1'  v  n  n  f 

(116) 

If  the  quantity  in  traces  in  Eq.  (ll6)  is  nonsingular,  then  it  has 
a  unique  inverse,  and  the  weighting  K(tn)  is  determined  uniquely  by 

K(t0)-P(t„/t„.l)if(tn){M(tll)P(tllAn.1)lf(t0)  *5 (t„)j  \  (117) 

Both  M(tn)P(t  /tn  ^M^t  )  &nd  R(tn),  which  is  defined  in  Eq.  (lOh), 

are  positive  semidefinite  from  their  definition  as  the  expected  values 

of  the  covariance  of  vector-valued  quantities,  so,  if  either  of  these 

two  matrices  is  positive  definite,  the  whole  quantity  in  braces  in 

Eq.  (112)  will  be  positive  definite  and  therefore  nonsingular. 

The  matrix  M(t  )P(t  /t  n )r  (t  )  will  be  positive  definite  if  the 

rows  of  M(tn)  are  independent  and  if  the  covariance  P(t  /t  ^)  is 

positive  definite.  The  covariance  P(t  /t  ,  )  is  positive  definite 

n'  n-1 

unless  there  is  some  x^(t  )  that  is  known  exactly;  if  it  were  known 
exactly;  it  wouldn't  have  to  be  estimated,  and  the  estimation  problem 
could  be  reduced  by  one  dimension.  When  the  rows  of  the  output  matrix 
M(t  )  are  dependent  and  there  is  no  multiplicative  noise  on  the 
measurements,  there  is  no  information  lost  by  reducing  the  number  of 
observed  random  variables  until  the  rows  of  M(t^)  become  independent. 
If  there  is  independent  multiplicative  noise  on  all  the  measurements, 
the  covariance  R(t^)  is  positive  definite  and  the  quantity  in  braces 
in  Eq  (ll6)  will  be  positive  definite.  Therefore,  even  if  the  quantity 
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in  braces  in  Eq.  (ll6)  is  not  nonsingular  to  start  with,  it  should  be 
possible  to  restate  the  estimation  problem  so  that  the  quantity  will  be 
nonsingular  in  the  new  statement  of  the  problem. 

If  for  some  reason  the  weighting  K(tr )  must  be  chosen  when  the 
quantity  in  braces  in  Eq.  (ll6)  is  a  singular  matrix,  the  weighting 
K(tn),  that  satisfies  Eq.  (ll6)  and  has  the  minimum  norm  of  all  the 
weighting  functions  that  satisfy  Eq.  (ll6),  is  given  by 

K(tn)  -  ♦  set, A  . 

(118) 

where  the  matrix  is  defined  as  the  psuedoinverse  of  the  matrix  A. 

A  particularly  concise  explanation  of  the  properties  of  the  psuedoinverse 
is  given  by  Gunckei  [Ref.  4]  in  an  appendix,  but  more  complete  explana¬ 
tions  can  be  found  elsewhere  [Ref.  26].  For  this  investigation,  it  is 
sufficient  to  say  that  the  psuedoinverse  gives  the  solution  to  a  set  of 
underspecified  equations  which  has  the  minimum  norm  of  all  possible 
solutions  to  the  set  of  equations.  In  the  remainder  of  this  paper  it 
will  be  assumed  that  the  matrix  in  braces  in  Eq.  (ll6)  is  nonsingular 
and  has  a  unique  inverse .  If  this  matrix  is  singular,  it  is  necessary 
only  to  replace  the  expression  for  the  inverse  of  this  matrix  by  the 
psuedoinverse  of  the  matrix,  and  all  the  equations  will  still  be  valid 
and  the  solution  will  have  the  minimum  norm  of  all  possible  solutions. 

E.  SUMMARY  OF  RESULTS 


The  procedure  developed  in  this  chapter  yields  an  iterative  solution 
to  the  problem  of  the  optimum  filtering  and  prediction  of  a  sampled  sto¬ 
chastic  process  that  can  be  represented  by  the  model  developed  in 
Chapter  II.  The  optimum  estimate  of  the  current  value  of  the  state 
vector  x('tn)  given  the  observed  random  variables  y(tn),y(tn  ^),..., 
y(tQ)  is 


x(t  /t  )  = 
x  n'  n 


I  -  ^n^nH  ^VVl^WVl)  + 


(119) 
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where  x(tn  ^/t  is  the  previous  optimum  estimate  and  the  ejected 

values  M(t  )  and  $(t  ,t  , )  are  known  a  priori.  The  r  x  p 
n'  v  n’  n-1'  — - - 

weighting  matrix  K(t^)  is  given  by 

K<‘»>  ■  +  ’ 

(120) 

where  the  covariance 

R(tn)  =  E[M(tn)x(tn)xT(tn)Mr(tn)]  (121) 


is  known  a  priori. 

The  covariance  matrix  of  the  error  in  prediction  P(t  ./t  )  is 
calculated  hy  the  recursion  relation 

P(t  ,/t  )  =  4(t  .  ,t  )[I  -  K(t  )M(t  )]P(t  /t  .  )^(t  ,  ,t  )  +  Q(t  ), 

v  n+1'  n'  v  n+1  n'L  v  n'  v  n/J  x  n'  n-1'  x  n+1'  n'  v  n" 

(122) 


where  the  covariance 


Q(tn)  =  E['.'>(tn+1,tn)x(tn)xT(tn)<I’T(tn+;L,tn)]  +  E[u(tn)uT(tn)]  (123) 

is  also  known  a  priori. 

Substituting  Eq.  (120)  into  Eq,  (122)  gives  the  nonlinear  difference 
equation  that  determines  the  covariance  of  the  error  as 


P(t  ,,t  ) 
v  n+1  n' 


It  is  interesting  to  compare  Eq.  (124)  with  the  equivalent  expres¬ 
sion  for  the  covariance  of  the  error  in  the  optimum  prediction  when  the 
parameters  of  the  process  are  not  random.  When  the  parameters  are  not 
random, 
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(125) 


M(t  )  -»  M(t  ) 

n  '  n' 

M(t  )  -  0 

n 

^Vl’V  -  <5(tn+l>tn) 

-*  ° 

R(t  )  -  0 

n 

Q(tn)  -*  Q(tn)  =  E[u(tn)uT(tn)], 

so  that  the  expression  equivalent  to  Eq.  (124)  is 

^  M(tn)j  P(t„/tn.1)*T(t11<.1,tn)  *  Q(tn).  (126) 

Equation  (126)  is  the  same  as  the  nonlinear  difference  equation  derived 
by  Kalman  [Ref.  1]  for  the  estimation  problem  when  the  parameters  oi  the 
process  are  known  a  priori,  instead  of  random  parameters.  The  effect 
of  the  random  variations  in  the  parameters  on  the  covariance  of  the 
error  in  the  optimum  prediction  P(tn  /t  )  ^-n  Eq.  (124)  is  included 

in  the  covariance  matrices  R(t  )  and  Q(t  ). 

When  the  parameters  of  the  process  M  and  i>  are  a  priori  known 
constants,  Kalman  [Ref.  l]  says  it  can  be  shown  that  the  optimum  estimate 
in  Eq.  (119)  becomes  a  stationary  dynamic  system  in  the  limit  as  the 
number  of  observed  random  variables  n  approach  infinity.  When  the 
matrices  M  and  and  the  covariances  R  and  Q  are  a  priori  known 
constants,  the  optimum  estimate  in  Eq.  (119)  should  also  become  a 
stationary  dynamic  system  in  the  limit.  The  reasoning  behind  the  latter 
statement  is  as  follows : 

The  trace  of  the  covariance  matrix  of  the  error  in  Eq,  (124)  is 
bounded  from  above  by  E[x(tn )xT(tn ) ]  (which  would  result  from 

an  estimate  that  the  state  vector  is  identically  zero)  and  from 
below  by  zero.  The  trace  for  the  optimum  estimate  should  not 
increase  from  one  estimate  to  the  next,  because  if  identical 
weightings  are  used  for  consecutive  estimates,  the  trace  of  the 
covariance  matrix  of  the  error  will  be  the  same  for  stationary 
statistics . 
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Thus,  the  trace  is  a  bounded  nonincreasing  sequence,  so  it  must 
approach  a  limit  as  n  approaches  infinity.  This  same  argument  can  be 
used  for  the  trace  of  the  covariance  of  the  error  in  the  optimum  esti¬ 
mate  of  any  linear  combination  of  state  variables;  therefore,  the  indi¬ 
vidual  elements  of  the  covariance  matrix  must  all  approach  a  limiting 
value.  The  limiting  value  of  the  weighting  matrix  K('tn)  is  obtained 
from  Eq.  (120),  therefore,  the  optimum  estimate  in  Eq.  (119)  becomes  a 
stationary  dynamic  system. 

The  matrix  block  diagram  in  Fig.  5  represents  the  linear  filter 

that  implements  the  optimum  estimate  in  Eq.  (119).  The  element  marked 

DELAY  relates  the  state  of  the  filter  at  time  t  .  to  the  state  of  the 

n+1 

filter  at  time  t  .  The  covariance  of  the  error  in  the  optimum  current 
estimate  in  Eq.  (119)  is  given  by 

P(tn/tn)  =  ^  '  K(tn)“(tn)1  P(tn/tn-l)  (l2?) 


so  that  from  Eq.  (124), 

OPTIMUM 

PREDICTION 


FIG.  5.  MODEL  OF  THE  OPTIMUM  FILTER. 


-  45  - 


P(t  ./t  )  =  0(t  ,,t  )P(t  /t  )^(t  ,,t  )  +  Q(t  ). 

n+1  n7  n+1  n7  '  n7  n7  v  n+17  n7  v  n7 


(128) 


The  optimum  prediction  of  the  state  vector  is 

x(t  +a/t  )  =  *(t  +a,t  )x(t  /t  )  a  >  o  (129) 

n  '  n7  v  n  ’  n7  v  n'  n7  —  v  7 

This  optimum  prediction  can  also  be  used  to  interpolate  between  sample 
points  when  the  a  is  time  varying  so  that 

a=t-t  t  <  t  <  t  , .  (130) 

n  n  -  -  n+1  v  J  7 

* 

The  initial  condition  of  the  linear  filter  with  no  observed  random 
variables  and  no  other  given  information  is  the  expected  value  of  the 
state  vector,  which  is  zero,  i.e., 

x(to/o)  =  0.  ;  (131) 

The  covariance  of  the  error  in  this  initial  estimate  is 

P(tQ/0)  =  E[x(tQ)xT(to)].  (132) 

The  weighting  K(tn)  is  determined  iteratively  at  each  sample  point 
from  Eqs ,  (120 )  and  (12*+)  and  the  initial  conditions. 

The  examples  in  the  following  section  will  illustrate  the  ideas 
of  this  chapter. 


F.  EXAMPLES 


The  procedure  for  estimation  and  prediction  derived  in  this  chapter 
uses  the  model  of  the  stochastic  process  developed  in  Chapter  II,  where 
the  stochastic  process  and  the  random  parameters  can  have  nonstationary 
statistics.  The  purpose  of  this  section  is  to  illustrate  and  clarify 
this  procedure;  therefore,  the  examples  will  be  limited  to  stochastic 
processes  and  random  parameters  with  stationary  statistics  For  the 
second-order  process  with  exponential  cosine  autocorrelation,  a  digital 
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computer  was  used  in  the  iteration  of  the  equations  which  determine  the 
weighting  function  and  the  covariance  matrix  of  the  error  at  teach  point. 

For  a  stationary,  stochastic  process  the  derivatives  of  the  covar¬ 
iance  of  the  state  vector  is  zero.  This  fact  leads  to  the  relation 
between  the  covariance  of  the  state  vector  E[x(t)x  (t)]  and  the 
covariance  of  the  input  V  where  F  and  G  are  matrices  defined  in 
Chapter  II. 

F  •  E[x(t)xT(x)l  +  E[x(t)xT(t)]FT  +  GVGT  =  0.  (133) 

The  covariance  Q(t^ )  can  be  determined  from  the  stationary  version  of 
Eq.  (112),  which  is 

E[x(tn)xT(tn)l  =  *(Tn)E[x(tn)xT(tn)]#(Tn)  +  Q(tn).  (134) 


The  solution  to  Eq.  (13*0  is  given  by 

r  r 


J=1  k=l 


1.  Exponential  Autocorrelation 

The  first  example  concerns  the  stationary  first  order  Markov 
process  which  was  presented  in  Fig.  2  and  discussed  in  Sec.  H-B.  The 
autocorrelation  function  of  the  process  is 

^(t)  =  e"plTL  (136) 

In  this  example  the  subscript  1  can  be  omitted  because  all  the  vectors 
and  all  the  matrices  are  scalars .  The  linear  differential  equation 
describing  the  system  is 


-  px(t)  +  V(t). 


(137) 


Therefore,  the  transition  matrix  is  a  scalar  given  by 
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(136) 


*  (Tn)  =  0(*n)  =  e  ^ 

The  expected  value  of  the  transition  matrix  ^(T^)  will  be  abbreviated 
so  that 

£  =  E[0{Tn)l  =  E[exp(-f3Tn)].  (139) 

From  the  autocorrelation  function  in  Eq.  (136)  the  covariance 
of  the  state  vector  E[x(t)x(t)]  is  given  by 


E[x(t)x(t)]  =  5^(0)  =  1. 


(140) 


Substituting  Eqs.  (139)  and  (l4o)  into  Eq.  (134)  gives  the  covariance 
Q(tn)  as 


Q(tn)  =  q(tn)  =  1  -  f-. 


(iki) 


The  multiplicative  noise  has  a  mean  of  unity  and  a  variance 
of  r,  so  that 


M(tn)  =  E[m(tn)l  =  1 

R(tn)  =  r  =  E[m(tn)x(tn)x(tn)m(tn)]  (142) 

The  model  of  the  optimum  filter  for  this  example  is  presented  in 
Fig.  6.  The  weighting  k(tQ)  Is  a  scalar  given  by  Eq.  (120)  as 


p(tn/tn  l') 

k(t  )  =  - 2 — Eli - 

p(t  /t  )  +  r 
'  n'  n-1 


(1^3) 


with  the  covariance  of  the  error  in  prediction  p(tn+^/t  )  calculated 
iteratively  from  the  nonlinear  difference  Eq.  (124),  which  is 


p(t  /t  ) 
n+1'  n 


0 


1  - 


P(tn/+n-l} 

^VVl)  +  ? 


p(tn/tn-l^  +  M2)*  (l44) 
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FIG.  6.  MODEL  OF  OPTIMUM  FILTER  FOR  A  STATIONARY  MARKOV  PROCESS. 

The  covariance  of  the  error  in  the  current  estimate  is  given  by 
Eq.  (127),  as 

p(t  /t  Jr 

P(t/t  )  =  - ;  n~L  _  •  (1^5) 

11  n  p(t  /t  J  +  T 

'  n'  n-1 

For  a  nonstationary  Markov  process  with  the  parameter  0  and  the 
covariances  ^  and  r  time-varying  quantities,  Eqs.  (ll+3),  (ll+l+),  and 
(ll+5)  would  still  be  valid,  but  0,  0,  and  r  would  have  different 
values  at  each  sample  point.  For  a  stationary  process,  the  covariance 
of  the  error  p(tn^/t^ )  will  approach  a  steady-state  value  as  n 
approaches  infinity. 

llm  P(tn+l/tn)  =  lirn  ^Wl)  =  P^n+l^n^  (lU6) 


The  steady-state  version  of  Eq.  (lUU)  is  a  quadratic  in  P(’tn+q/t,n)> 


-  1+9  - 


I 


p(t  .,/t  )  = 
'  n+17  n 7 


$2  r  p(t  /t  ) 

- n+l_n_  +  (1  _  f 2 

p(t  ../t  )  +  r 
v  n+17  n7 


(l^T) 


which  h..s  the  solution 


A-A’  ■  I  h  -  f) 


(l-r)  +,/  (1-r)  + 


— x2  uT"1 


(1U8) 

(The  square  root  in  Eq.  (l48)  has  a  positive  sign  because  p(t  /t  ) 
must  always  be  positive, ) 

Three  common  proabability  laws  for  the  sample  period  T^  were 
discussed  in  Sec,  II-E.  The  pertinent  information  for  these  three 
laws  for  particular  numerical  values  of  the  parameters  of  the  proba¬ 
bility  laws  and  for  periodic  sampling  is  presented  in  Table  2.  The 
expected  value  of  the  sample  period  is  normalized  to  unity,  and  the 
variances  of  the  sample  period  have  the  values  stated.  The  expected 
value  of  the  transition  matrix  0  is  listed  for 


P  =  0.1. 


(149) 


For  all  numerical  calculations  this  value  of  (3  will  be  used;  therefore, 
the  system  will  have  a  time  constant  of  ten  seconds.  Notice  how 
increasing  the  variance  of  the  sample  period  increases  the  expected 
value  of  the  transition  matrix  0. 

The  steady-state  values  of  the  covariance  of  the  error  P(tnAn) 
and  p(t  ^/t  )  are  presented  in  Table  3  for  the  same  expected  values 
of  0  listed  in  Table  2  for  two  cases — with  no  multiplicative  noise  and 
with  the  variance  of  the  multiplicative  noise  equal  to  one. 

Notice  that  with  optimum  estimation  the  variation  in  the  proba¬ 
bility  distribution  of  the  sample  period  has  a  smaller  effect  on  the 
covariance  of  the  error  than  does  the  variance  of  the  multiplicative 
noise  The  two  covariances  PAnAn)  and  p(t  ^/t  )  are  plotted  in 
Fig.  7  as  a  function  of  the  variance  of  the  multiplicative  noise  r 
for 
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TABLE  2.  NUMERICAL  VALUES  FOR  FOUR  PROBABILITY  LA*S. 


Purely  Randoa  or  Poieion  Stapling 


e  Sampling  with  Independent  Miaaea 


Periodic  Sampling  with  Jitter 


Periodic  Sampling 


probability 
Di  a t r i bu ti on  of 
Sample  period 


Parameter#  of 
pro  babi 1 i ty 
Distribution 


Mean  Sample 
Period 

E[Tn] 


Variance  of 


Samp 1 e  period 
E[T  2]  -  (  E  [T  ]  )  2 


4>  ”  E  [exp(  )  ] 
10 


TABLE  3.  COVARIANCE  OF  THE  ERROR  FOR  EXAMPLE  1. 


pro  babi 1 i ty 
Distribution 

o  f 

Samp  1 e  period 

Co n  s  t  an  t 

Normal 

Geom  e trie 

Expon  en  ti al 

r  T 

<£-E  [exp(  --ft  ) 
lu 

0.9048 

0.9053 

0. 9070 

0.9091 

p(tn/tn^  with 

r  =0 

0.  0 

0.0 

0  .  0 

0.  0 

p(t  + 1  / 1  )  with 

n  n 

r  =  0 

0.  18  13 

0.  1804 

0.  1773 

0. 1736 

P(tn/tn)  with 

r  =  1 

0. 2986 

0.2982 

1 

0.2962 

0. 2927 

P<  tn+1/tn> wi th 

r=  1 

0. 4258 

! 

j  0.4247 

0.4211 

0.4167 

L 
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FIG.  7.  COVARIANCE  OF  THE  ERROR  FOR  EXAMPLE  1. 

?  =  0.9091,  (150) 

whicn  is  the  case  where  the  sample  period  has  an  exponential  probability 
distribution  (purely  random  sampling). 

2.  Exponential  Cosine  Autocorrelation 

The  second  example  concerns  the  stationary  process  with  the 
exponential  cosine  autocorrelation 


31 


X1X1 


(t) 


P 7  sin  7 
2P2  +  7 


(151) 


that  was.  presented  in  Fig.  3  and  discussed  in  Sec.  II-C.  The  linear 
differential  vector  equation  that  describes  the  system  is 


d_ 

Xx(t) 

-  p  1 

x1(t) 

+ 

vL(t) 

dt 

Xg(t) 

CQ_ 

1 

OJ 

1 

x2(t) 

0 

—  — 

L 

_ 

- 

(152) 
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The  transition  matrix  is  given  by  Eq.  (40)  as 


cos  y  T  —  sin 

'  n  •  y  n 


-y  sin  -vT  cos  >T 
'  '  n  '  n 


When  the  sample  period  has  a  geometric  distribution,  the  expected  value 
of  the  transition  matrix  can  be  calculated  by  changing  the  sines  and 
cosines  to  complex  exponentials,  so  that 


-pT 

E[e  cos  ?T  ]  =  RE [exp  (-PTn  +  JyT  )]  =  2 

(d  +  P)  +7 


E[e  n  sin  yT  ]  =  lM[exp  (-PT  -  J7T  )]  =  - - 2  (15*0 

(d  +  P)  +  7 


where  RE  and  1M  represent  the  real  and  imaginary  parts,  respectively. 

T 

The  covariance  of  the  state  vector  E[x(t)x  (t)]  is  determined 
from  Eq.  (133).  For  this  example  the  solution  to  Eq.  (133)  is 


E[x1(t)x1(t)]  =  Sx  x  (0)  =  1 


E[xL(t)x2(t  )]  =  E[x2(t)x1(t)]  =  -  p  y2/{y2  +  2p2) 


E[x2(t)x2(t)]  =  yk/{y2  +  2p2) 


E[v1(t)v1(t)]  =  VL1  =  4  p  (y2  +  p2)/(72  +  2P2).  (155) 

The  covariance  Q('tn)  is  determined  directly  from  Eq.  (135)-  The 
multiplicative  noise  has  a  mean  of  unity  and  a  variance  of  r,  so  that 

E[“T(t»h  ■  [0] 

R(tn)  =  E[M(tn)x(tn)xT(tn)Mr(tn)]  =  r  E[x1(tn)x1(tn)]  =  r.  (156) 
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The  numerical  values  of  the  transition  matrix  o(T  )  and  the 
covariances  Q^t^)  and  E[x(t)x  (t)]  are  presented  in  Table  4  for  a 
geometric  distribution  of  the  sample  period  (purely  random  sampling) 
with  y=l,  so  that 


H  =  1.0 

P  =  0.1 

y  =  1.0. 


(157) 


TABLE  4.  NUMERICAL  VALUES  OF  MATRICES  FOR  EXAMPLE  2. 


$(TJ 

0.4977  0.4525 

l_-  0  .  4  52  5  0.  4977 

6u„) 

0.5956  -0.08964 

-0.08964  0.4886 

E[x(t)xT(0] 

1.0  -0.09804 

-0.09804  0.9804 

The  expected  value  of  the  transition  matrix  $(T  )  and  the  covariance 
_  n 

Q(t^)  are  needed  for  the  iterative  calculation  of  the  covariance  of  the 
error  P(tn  /t  )  from  the  nonlinear  difference  Eq.  (124).  The  initial 
condition  for  the  covariance  of  the  error  P(tQ/0)  with  no  other 
information  is 

P(tQ/0)  =  E[x(t)xT(t)].  (158) 

Because  the  observed  random  variable  is  a  scalar,  the  quantity  in  Eq. 
(124),  which  must  be  inverted,  is  also  a  scalar. 
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l“  (V^Wi^V  *  A*)'1  ■  vtWi/V  * F1  (159) 

For  this  example  the  steady -state  covariance  of  the  error 
P(t  ^/tjj)  was  obtained  by  iteration  of  Eq.  (124)  on  a  digital  computer. 
The  steady -state  value  of  the  covariance  P(t  ^/t^)  and  the  weighting 
matrix  K(tn)  are  presented  for  the  numerical  values  of  the  matrices 
in  Table  5  when  the  multiplicative  noise  has  a  variance  of  one -tenth, 

r  =  0.1.  (160; 

TABLE  5.  DESIGN  RESULTS  FOR  EXAMPLE  2 


P^n+1/t„) 


K(t„) 


0.  7  56  4 

0.0411 

0  0411 

0.6700 

K(t0) 


K  (  t  j  ) 


K(t2) 


K(  t3) 


0.  88  32 
0.  0480 


0.909  1 
0.0893 


0.8905 
0  1190 


0.  8  8  56 
0  .  06  56 


0  8  8  37 
0  05121 


The  weighting  matrix  K(t  )  is  also  presented  j.n  Table  5  for  the  first 
four  sample  points  to  show  hew  fast  the  iterations  converge  to  a  steady - 
state  value . 

From  Eq.  (127)  the  steady-state  value  of  the  covariance  p^  (tn/tn) 
is 


pU(tn+A>  r 


P..(t  /t  )  -  " - =.  (163 

p  (t  /t  )  +  r 
11 v  n+1'  n 

In  Fig.  8  the  steady-state  covariances  of  the  error  p...(t  /t  )  and 


0.7 


VAR  i  AHCE  Of  THE  NODE  7 


FIG  8  COVARIANCE  OF  THE  ERROR  FOR  EXAMPLE  2. 

P1]L(tn+^/t  )  are  plotted  as  a  function  of  the  variance  of  the  multi¬ 
plicative  noise  r  for  a  geometric  distribution  of  the  sample  period 
for  the  value  of  the  matrices  in  Table  4.  Notice  how  much  the  cosine 
term  in  the  autocorrelation  has  increased  the  covariance  of  the  error 
in  Fig.  8  over  the  corresponding  covariances  in  Fig.  7- 
3.  Nonstationary  Process 

The  third  example  concerns  a  first-order  Markov  process  with  a 
particular  kind  of  nonstationarity.  The  linear  differential  equation 
describing  the  system  is 

^  =  -  P(t)x(t)  +  v(t) 


with  p(t)  time  varying  so  that 


p(t)  =  3  2k  <  t  <  2k  +  1 


=B2  2k  +  1  <  t  <  2k  k  -  0,1,2,  . . 
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Given  the  value  of  the  process  at  time  2k  +  t  with 


0  <  T  <  2, 

it  is  desired  to  predict  the  value  of  the  process  at  time  2k  +  T  +1. 
From  Eq.  (l29)  the  prediction  is  obtained  by  multiplying  the  current 
value  of  the  process  x(2k  +  T )  by  the  transition  matrix  0(2k+T+l,2k+r) 

x(:2k+T+  l/2k+T)  =  0(2k+T+  '1, 2k+T  )x(2k+T ) 

-B,T  -B2T2 

=  e  1  1  e  2  2  x(2k+T ) 


where 

T  =  1-T  Tg  =  t  for  0  <  T  <  1 

T  =T  -1  T2  =  2-t  for  1  <  T  <  2 

If  v(t)  has  unity  variance,  the  mean  squared  error  of  the  estimate 
is 


-2B2T2  1  -  e 


-2^1 


1  -  e 


-2B2T2 


2B, 


2B. 


These  examples  have  showed  problems  in  prediction  of  random  processes 
„hat  can  be  solved  by  using  the  method  developed  in  Chapter  III. 
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IV.  OPTIMUM  INTERPOLATION  WITH  DELAY  OF  RANDOM  PARAMETER  PROCESSES 


A.  INTRODUCTION 

In  Chapter  III  the  matrix-valued  sampledWiener-Hopf  equation  was 
derived,  and  the  optimum  estimate  of  the  current  value  of  the  state 
vector  vas  shovn  to  be  a  linear  combination  of  the  previous  optimum 
estimate  and  the  nevly  observed  random  variable  at  the  current  sample 
point.  In  this  chapter  it  vill  be  proved  that  any  optimum  estimate  that 
satisfies  the  sampled  Wiener-Eopf  equation  can  alvays  be  vritten  in  an 
iterative  form  similar  to  the  one  derived  in  Chapter  III.  The  exact 
configuration  of  this  iterative  form  is  worked  out  for  the  general  pro¬ 
blem  of  interpolation  with  delay  (at  time  t^  the  optimum  estimate 

is  desired  of  the  state  vector  at  time  t  where  t  <  t  ).  and  the 

n  7 

iterative  nonlinear  difference  equations  are  derived  for  the  matrix¬ 
valued  gain  of  the  optimum  filter  with  delay  up  to  (n  -  d)  sample 
periods.  The  reason  for  delaying  the  estimate  is  that  the  trace  of 
the  covariance  matrix  of  the  error  of  the  optimum  estimate  can  be 
decreased  when  the  estimate  is  delayed  until  additional  random  variables 
are  observed. 

B.  ITERATIVE  SOLUTION  TO  THE  SAMPLED  WIENER-HOPF  EQUATION 


In  this  section  it  will  be  proved  that  at  any  sample  point  t 
the  optimum  estimate  x(t  +a/t  )  can  be  written  in  the  iterative  form 


x(t  +a/t  )  =  x(t  +a/t  n)  +  K(t  +a,t  )[y(t  )-M(t  )x(t  /t  _,)] 
'  n  '  n  '  n  n-1  n  n  L  n  '  n  n'  n-l'J 


(162) 

where  x(t  +a/t  is  the  optimum  estimate  at  the  previous  sample  point 

and  K(t  +a,t  )  is  an  r  x  p  weighting  matrix.  When  a  is  positive, 
n  n 

the  estimate  x(t  -Kj/t^)  is  called  prediction,  and  when  a  is  negative, 
it  is  interpolation  with  delay.  The  proof  in  this  section  is  a  general¬ 
ization  of  the  proof  in  Chapter  III  that,  when  a  is  zero,  the  optimum 
current  estimate  can  be  written  in  the  form  of  Eq.  (162) , 
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Because  the  two  estimates  x(t  +a/t  , )  and  x(t  /t  ..  )  are 

'  n  '  n-1'  n'  n-1' 

optimum,  they  satisfy  the  appropriate  versions  of  the  sampled-Wiener- 
Hopf  equation  which  are 

E[x(tn+0!)yT(tp)]-E[x(tn+a/tn_1)yT(t  =0  p  =  0,1,  . ..,n-l 

(163) 

and 

E[x(tn)yT(tp)]-E[x(tn/tn_1)y'I'(tp)]  =0  p  =  0,1, . . .  ,n-l. 


(164) 

It  will  be  proved  that  there  is  an  r  x.  p  weighting  matrix  K(t^+a,tn) 
for  which  the  estimate  in  Eq.  (l62)  satisfies  Eq.  (165), 


E[x(tn+o:)yT(t  )]-E[x(tn+a/tn)yT(t  )]  =0  p  =  0,l,...,n. 

(165) 

Substituting  the  trial  estimate  in  Eq.  (162)  into  Eq.  (165)  gives 


Erx(tn+a)yT(tp)l-E[x(tn^/tn_1)yT(tp)l 


-K(tn+a^tn)  ^E[M(tn)x(tn)yT(tp)l 
-M(tn)E[x(tn/tn_1)yT(tp)]^  =  0 


p  —  0^1^ 


(166) 

For  p  less  than  n,  Eq.  (l66)  is  equal  to  Eq  (163)  minus 
K(tn+a,tn)M(tn)  times  Eq.  (l64).  Both  these  equations  are  identically 
zero  for  p  <  n,  so  Eq.  (l66)  is  identically  zero  for  p  <  n. 

For  p  equal  to  n,  Eq.  (l66)  can  be  written 


E[x(tn+a)xT(tn)]-E[x(tn+a/tn_1)xT(tn)] 


f?(tn) 


—I 

-K(t  +a,t  )  l  M(t  )p(t  /t  .Jlf^t  )+R(t  )\ 

n  n  1  n'  n'  n-1'  '  n  n  | 


=  0, 


(167) 
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where  the  substitution  has  been  made  that 


E[M(tn)xl(trl)yT(tri)T-E[x(tn/tn_1)yT(tn) 

=  M(tn)P(tri/tn_1)Mr(tn)+R(tri),  (168) 

with  the  covariance  matrix  of  the  error  P(t  /t  ,  )  defined  as 

n'  n-1 

P(tn/tn-l^  =  E[x(tn)xT(tn)1-E[S(tn/tn.1)xT(tr_)]>  (169) 

and  the  covariance  matrix  R(tn)  defined  in  Eq.  (104)  as 

R(tn)  =  E[M(tri)x(tn)xT(tn)lf  (tn)],  (170) 

Therefore,  when  the  rxp  weighting  matrix  K(tn+a,tn)  satisfies 

Eq.  (167),  the  trial  estimate  in  Eq.  (162)  is  indeed  optimum  and  satis¬ 
fies  the  sampled-Wiener-Hopf  equation  (165).  When  the  quantity  in 
braces  in  Eq  (167)  is  nonsingular,  the  weighting  matrix  K(tn+a,tn) 
is  determined  uniquely  as 

K(t  +a,t  )  =  <E[x(t  +a)xT(t  )]-E[x(t  +a/t  )xT(t  )]  \  i?(t  )  • 
v  n  n  |  1  n  '  n  '  1  n  '  n-1  n  J  I  '  n 

{S(tn)P(tn/tn-l)”I'(tn)+I(tn)}  0-71) 

The  conditions  under  which  the  second  quantity  in  braces  is  nonsingular 
are  discussed  in  Sec.  III-C.  When  a  is  equal  to  zero,  the  weighting 
K(tn,tp)  ^or  optimum  current  estimate  from  Eq.  (171)  is 

K^tn,tn^  =  ^E[x(tn)xT(tn^_E[5(tn/tn-l^xT^tn^j  ^0^)  ’ 

^M(tn)P(tnAn_i)Mr(tn)+R(tn)j  (172) 

which  is  the  same  as  K('tn)  in  Eq.  (120). 

I 
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It  is  also  desirable  to  know  the  covariance  matrix  of  the  error  of 
the  optimum  estimate  P(tn+a/tn),  which  is  given  by  Eq.  (87)  as 


P(tn+0[/tn)  =  Etx(tn4°!)xT(tn'H:it)]  '  E[x(tnWtn)xT(tn+0!)). 


(173) 


Substituting  the  optimum  estimate  in  Eq.  (162)  into  Eq.  (l7l)  gives 

p(tn+a/tn)  =  E[x'.t„+a)xT(tr,+a)l  -  E[x(t„+a/tri_1  )xT(tr(+a)] 


n-ly 


-  K(tn+a,tn)M(tn)^E[x(tn)xT(tn4a)]-E[x(tn/tn_1)xT(tn+a)]j 

(175) 

The  first  two  terms  in  the  expression  for  the  covariance  P(t  +a/t  ) 
in  Eq.  (175)  are  equal  to  the  previous  covariance  P(t^+a/t  ^);  thus, 
the  covariance  matrix  of  the  error  of  the  optimum  estimate  can  also  be 
written  as 

P(tn+a/tn)  =  P(tnW'tn_1)-K(tn+a1tn)M(tn)  ^E[x(tn)xT(tn+a)] 

-  E[x(tn/tn.1)XT(tn+a)1^ 


(176) 


When  a  is  equal  to  zero,  the  covariance  Pft^-Hi/t^)  in  Eq.  (176)  is 
equal  to  the  covariance  p(t  /t  )  in  Eq ,  (127) 

In  the  next  section  the  iterative  nonlinear  difference  equations 
determining  the  weighting  matrix  K(t^+a,tn)  in  Eq  (172)  and  the 
covariance  matrix  P(t  -Kl/t^)  in  Eq.  (176)  will  be  derived  for  inter¬ 
polation  with  delay  (when  a  is  negative). 


C.  DERIVATION  OF  THE  OFTIMUM  ESTIMATE 


To  simplify  the  derivation  of  the  optimum  estimate  for  interpolation 
with  a  delay  of  |a|  sec,  the  notation  used  in  the  preceding  section 
will  be  changed  slightly  The  time  t^+a  can  be  written 
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(177) 


t  +  a  =  t ,  +6  0  <  £>  <  T , 

n  d.  —  a 

so  that,  for  a  negative, 

t  <  t  +  a  <  t .  ,  a  <  0 

*  —  n  d+1  — 

d  <  n.  (178) 

Using  the  new  notation  the  optimum  estimate  x(t  +a/t  )  in  Eq„  (162) 
becomes 

x(td+b/tn)  =  x^-ts/t^)  +  K(td4©,tn)[y(tn)-M(tn)x(tn/tn_1)].  (179) 

Making  use  of  Eq,  (l62)  once  more  the  optimum  estimate  x(td+&/t  ^ )  in 
Eq.  (179)  can  be  written 

^V^n-l*  =  &(td^/tn-2)+K(td^’tn-l)fy(tn-l)-M(tn-l)i(tn-l/tu~2)^ 

(!8°) 

Continuing  this  process  for  (n  -  (■)  times  the  expression  for  the 
optimum  estimate,  x(td+&/tn)  becomes 

n_^ 

x(td4b/tn)  =  xft^/t^)  +  ^  K(td+6,tk)[y(tk)-M(tk)x(tk/tk_1)].(l8l) 

k=i+l 

For  a  delay  of  (n-d)  sample  points,  it  is  necessary  to  calculate  (n-d) 
weighting  matrices  K(td+&,tk). 

The  matrix  block  diagram  in  Fig,  9  represents  the  linear  filter 

which  implements  the  optimum  estimate  in  Eq,  (l8l).  The  elements  marked 

DELAY  relate  the  state  of  the  filter  of  time  t  to  the  state  of  the 

n 

filter  at  time  t  That  part  of  the  filter  within  the  dashed  rec¬ 

tangle  duplicates  the  model  of  the  optimum  filter  for  prediction  in 
Fig.  5>  which  was  derived  in  Chapter  III,  and  calculates  the  optimum 
prediction  x(t  /t  )  in  Eq.  (l8l). 

K  K  ~ 

The  weighting  matrix  K(td+6,tn),  where  n  represents  any  number 
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FIG.  9.  MODEL  OF  OPTIMUM  FILTEH  WITH  DELAY. 

greater  than  d,  must  satisfy  Eq.  (167)  with  t^+a  replaced  by 
td+s;  i.e., 

|E[x(td+6)xT(tn)]  -  E[x(td46/tn_1)xT(tn)]^MT(tn) 

-  K(td+6,tri)[M(tn)P(tn/tn_1)Mr(tn)TR(tn)]  =  .0. 

(182) 

The  covariance  of  the  error  PA^A^-p)  in  E<1>  (l82)  can  be  determined 
iteratively  at  each  sample  point  by  Eq.  (124);  thus,  only  the  first 
term  in  braces  in  Eq.  (l82)  remains  unknown.  This  term  is  composed  of 
the  sum  of  two  quantities  E[x(td+6  )xT(tn) ]  and  E[x(td+&/tn  d)xT(tn)]. 
This  sum  will  be  determined  first  for  the  case  where  there  is  a  delay  of 
only  one  sample  point,  so  that  n  is  equal  to  d+1.  A  method  will  then 
be  derived  for  calculating  this  sum  in  the  general  case  where  there  is  a 
delay  of  (n-d)  sample  points. 

From  the  solution  to  the  ordinary  differential  equation  representing 
the  system  in  Eq.  (10),  the  state  vextor  x(td+6 )  can  be  written 
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x(t,46)  =  *(t.-t6,t,)x(t,)+u(t.+6,t,). 


(183) 


where  u(td+6,td)  is  an  r  x  1  vector  defined  as 

td+B 


u(td+6,td)  =  J  4(td4€,,T)G(T)v(T)dT 


(184) 


to  represent  the  effect  of  the  random  inputs  to  the  system  between 
time  t  and  time  td+6 .  From  the  model  of  the  process  and  the 
properties  of  the  random  parameters  in  Eq.  (53)  the  state  vectors  can 
be  written 


x(td46)  =  *(td+6,td)x(td)  +  u(td4S,td) 


=  ®(td4b,td)x(td)  +  $(td46,td)x(td)  +  u(td-tb,td)  (185) 


and 


x(td+i)  =  «(Vl>td)x(td)  +  u(td+l'td) 


=  ®(td+ljtd)x(td)  +J(td+l'td)x(td)+u(td+l'td)‘  (l86) 
Substituting  Eqs.  (l85)  and  (l86)  into  the  quantity  E[x(t  +fj)xT(t  ,  )) 


gives 


E[x(td4b)xT(td+1)]  =  $(td4$,td)E[x(td)xT(td)] 

(187) 

where  the  expected  values  of  the  cross  terms  are  zero,  and  the  substi¬ 
tution  has  been  made  of  the  r  x  r  matrix  Q^Ct^+g)  defined  as 

Qx(td4b)  =  E[5(td+5,td)x(td)xT(td)?T(td+1,td)]+E[u(t(i^,td)uT(td+1,td)], 


(188) 


-  6k  - 


The  quantity  E[x(td4$/td)xT(td+d)]  can  be  written 

E[x(td+&/td)xT(td+1)]  =  $(td-«,td)E[x(td/td)xT(td)]  ^r(td+1>tdL)‘  (1®9) 

Thus,  by  substituting  Eqs.  (187)  and  (189),  the  sum  of  the  two  quantities 
E[x(td+5)xT(td+1)]  and  E[x(td4$/td)xT(td+1)]  can  be  written 

E[x(td+5)xT(td+1)]  -  E[x(td+&/td)xT(td+1)] 

=  ^(td+&,td)E[x(td)xT(td)]fT(td+1,td)-^1(td+6) 

-  ^td+&'td)E[S(td/td)xT(td)]5T(td+l'td) 

=  ®(td+&,td)p(td/td)  ®r(td+1,td)+Q1(td4e),  (190) 

where  the  substitution  has  been  made  that 


P(td/td)  = 


d'  d' 


(191) 


Returning  to  the  general  case  with  the  first  term  in  braces  in 

T 

Eq.  (182)  composed  of  the  sum  of  E[x(td+&)x  (tn)]  and  E[x(td+fj/tn  d) 
xT(tn)],  the  substitution  of  Eq.  (l8l),  for  x(td+&/tn  gives 


E[x(td+©)x  (tn)]  -  E[x(td4f,/t. 


=E[x(td+6)xT(tn)]  -  E[x(td4h/td)xT(tn)] 
n-1 


'  £  (E[y(tk)xT(tn.i)]-M(tk)E[x(tk/tk-l)xT(tn)l'] 

k=d+l  ^ 


For  k  less  than  n, 


(192) 


{E[y(tk)x'r(tn_1)]  -  M(tk)E[x(tk/tk_1)xT(tn)]] 

^  _  r  m  m 


=  M(tk)J  EjXt^x1^)]  -  E[x(tk/tk_1)xJ’(tk)]  |  ^T(tnJtk) 


-  ►"‘ki'VVi' 


(193) 
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where  the  covariance  of  the  error  P(t^/t^  )  has  been  substituted  into 

Eq.  (193)*  For  d  less  than  n 

E[x(td  +&)xT(tn)]  -  E[x(td+&/td)xT(tn)] 

*  (♦<td+6'td)p(td/td)«T(tati>td)+\<ta,6)}  ?(VVi)  (W1*) 

where  Eq.  (190)  has  been  substituted  into  Eq.  (19^). 

Finally,  substituting  Eqs.  (192  -  19^)  into  Eq.  (l82),  which  must 
be  satisfied  by  the  optimum  weighting  matrices,  yields 

ptd**"td>P<td/td)#(td.l’td)+5l<ta,6))*T<td>td+l)i?<tn> 

-  X  K(td«,tlt)»(tk)P(tk/tk_1) 

■  k=d+l 

’  K(td+&>tn)  =  °’  (195) 

where  )  and  R(t  )  are  defined  in  Eqs.  (l88)  and  (167).  First, 

set  n  equal  to  d+1,  and  Eq.  (195)  becomes 

r*<td*,.dW*dAd>  ^>j  if'<tn> 

-K(td4t,,tdti)r«(Vi)P(tdti/td)Mr(td+i)*B(td+i)j  -0.  (196) 

The  optimum  weighting  matrix  K(td+6,td+d)  Is  determined  from  Eq.  (196) 
to  be 

‘(V'Vl1  =  (*<‘d^d)P<V‘d)  *'W<5l(*d*)}if(t4) 

_1-  <W) 

Because  the  delay  5  is  known,  the  expected  value  o(td+6,td)  and 
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can  be  calculated  directly;  thus,  all  the  quantities  in  Eq.  (197) 
are  known  a  priori  except  for  the  covariances  of  the  error  P(t^+^/t^) 
and  P(t^/t^),  which  can  be  determined  iteratively  at  each  sample  point 
by  Eqs.  (124)  and  (127).  The  conditions  under  which  Eq.  (197)  uniquely 
determines  the  weighting  manrix  K(t^+&,td+1)  are  discussed  in 
Section  III-D. 

Next,  set  n  equal  to  d+2  in  Eq.  (195)  and  the  result  is 

j^<V6,ta)P(td/td)  ^(ta+1,ta)*5(t1«)j  jT(V?'td+l)i?'(td*2) 

-  lt<ta*,t1+1)M(ta+i)p(ta+iAa)  5T<V2’Vi>if<ta+2> 

-  K'td«.‘d+2>(R<td+2>p<ta+2''tdTi>i?(td+2>+*<V2>}  ■  °- 

(198) 

The  optimum  weighting  matrix  K(td+B>td+2)  is  determined  from  Eq.  (198), 
to  be 

*(t«<e’td+2)  ■  {!F<td«,ta)P(taAa)  ®r(ta+2,ta)Mr(ta+2) 

A<V*> 

-K(ta+B,td+i)if(td+1)P(ta+1Aa) 
{R(ta.2)P(V2/tdti>if<V2>+E<td+2)j  _1  (199) 

with  K(td+6»td+1)  given  by  Eq.  (198). 

To  calculate  the  optimum  estimate  x(t,+6/t  )  in  Eq.  (l8l)  with 
delay  up  to  (n-d)  sample  periods,  it  is  necessary  to  determine  (n-d) 
weighting  matrices  K(t,+6,t  ).  The  weighting  matrices  K(t  +6,t  ) 

Cl  K.  Q.  CItJ. 

and  K(t^+6,t^+;j)  have  already  been  determined  from  Eq.  (195)  as 
given  by  Eqs.  (197)  and  (199)  respectively.  The  matrix  K(t^+&,t^+^) 
is  obtained  from  Eq.  (195)  with  n  equal  to  3>  and  this  iterative 
procedure  is  continued  until  all  the  (n-d)  weighting  matrices  have 
been  obtained 
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These  calculations  will  be  rather  tedious  if  the  optimum  estimate 
has  been  delayed  fc1”  many  sample  points  and  there  are  a  large  number 
of  weighting  matrices  to  be  determined.  When  the  optimum  estimate  is 
obtained  from  a  stationary  dynamic  system  (as  discussed  in  Sec.  III-E), 
these  weighting  matrices  will  not  change  from  one  sample  point  to  the 
next;  thus,  the  tedious  calculations  must  be  performed  only  once. 

These  ideas  are  illustrated  in  Sec.  IV-D,  where  examples  are  worked 
out  for  the  optimum  estimate  with  a  delay  of  one  sample  period  and  with 
a  delay  of  two  sample  periods.  For  a  delay  of  up  to  one  sample  period 
the  covariance  of  the  error  P(t^+b/t^+^)  is  given  by  Eq.  (175),  with 
t^+a  and  t  replaced  by  t^+6  and  t  respectively,  so  that 

P(td+&/td+1)  =  E[x(td+6)xT(td+6)]-E[x(td+&/td)xT(td-«)] 

"K  ( 1  d+6 ,  Vi )R(  td+i }  (E  [ x  ( Vi  )xT(td+&)] 

-E[J(td+1/td)xT(td+6)]j  .  (200) 

T 

By  using  Eq.  (185)  the  covariance  E[x(td+b)x  (td+b)]  can  be 
written 

E[x(td+6)xT(td+b)l  =  ^(td+6,td)E[x(td)xT(td)]  ^I'(td+6,td)+Q2(td+6), 

(201) 

where  the  expected  value  of  the  cross  terms  is  zero,  and  the  substitution 
has  been  made  of  the  r  x  r  matrix  Q2(td+&)  defined  as 

Q2(td+&)  =  E[^(td4b,td)x(td)xT(td)?T(td+b,td)]+E[u(td-^,td)uT(td+6,td)] 

(202) 

Both  the  quantity  Q2(td+6)  in  Eq.  (202)  and  the  quantity 

defined  in  Eq.  (188)  are  equal  to  zero  when  t  +  5  is  equal  to  t  . 

d  d 

Substituting  Eqs.  (187)  and  (201)  into  Eq.  (l99)>  and  using  the 

covariance  P(td/td)  for  the  equivalent  quantity  given  by  Eq.  (191) 

yields  as  the  covariance  matrix  of  the  error  P(t  +&/t  ) 

a  d+1 
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P(td46/td+l)  =  ^(td46,td)P(td/td)  ^^td4*>td)  + 

-  K<td*«W>“<Vi>f' *<Vi',a>p(‘dAd) 

+  ^<ta«A  ,  (203) 

where  K(td+5,td+1)  is  given  by  Eq.  (197),  and  and  Q2(td+&) 

are  defined  by  Eqs.  (l88)  and  (202)  respectively. 

For  a  delay  of  up  to  two  sample  periods  the  expression  for  the 
covariance  of  the  error  P(td+g Ad+2)  can  be  derived  from  Eq.  (176) 
by  the  same  procedure  used  to  derive  P(td+8 Ad+d )  from  Eq.  (176). 

P(td46/td+2)  =  P(td+e/td+l)-K(td4*’td+2)^(td+2)(Etx(td+2)xT(td46^ 

-E[S(td+2/td+i)xT(td-«)]j  •  (204) 

Everything  in  Eq.  (204)  has  already  been  determined  except  for  the 
quantity  in  braces.  From  Eq.  (162)  the  optimum  estimate  x(t,  ,,/t  ,) 

G.+<^  a+j_ 

can  be  written 

x(td+2/td+l)  =  ^WW^WVl* 

=  ^td+2,td+l^x^td+l//td^+K^td+l^y^td+l^‘^  ^d+l^  x^td+l^td^  J 

(205) 

Substituting  Eq.  (205)  into  the  quantity  in  brackets  in  Eq.  (204)  gives 
j^E  [  X(td+2  )XT  ( td+B  )  ] -E  [  x(ti+2/td+i  )xT  ( td4f, ) 

=  -  «tdti«ViM  x 

|^E[x(td+1)xT(td+8)l-E[x(td+1/td)xT(td+&)]^  .  (206) 

Substituting  Eqs.  (206)  and  (187)  into  Eq.  (204)  gives  as  the  covariance 
of  the  error 
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P(V6/td+2>  ■  P(td*/td+l)-K(td«'td+2)M<td»2>  ‘<V2'W>X 


[I-K(Vi)H(td+i)l[  ♦(td+1,ta)P(t4/t4) 

(207) 

where  I  is  the  identity  matrix  and  P(t^46/t^+^)  and  Kit^-te^t  g) 

are  given  by  Eqs.  (203)  and  (199)  respectively. 

For  a  delay  of  up  to  (n-d)  sample  periods  the  general  expression 

for  the  covariance  matrix  of  the  error  is  given  by  Eq.  (176),  and,  by 

repeating  the  same  procedure  used  to  obtain  Eq.  (205),  this  covariance 

P(t ,+6 /t  )  can  be  written 
'  d  '  n 

n-1 

P(td+6/tn)  =  P(td+6/tn-l)+K(td+6^tn)  E  *(tk+l'V [' I"K(tk)M(tk) ]  * 

k=d+l 

[  ®(td+l’td)P(td/td)  (208) 

Therefore,  the  covariance  of  the  error  ^(t^-fs/t^)  can  be  calculated 
by  starting  with  Eq.  (203)  and  using  Eq.  (208)  iteratively. 

In  conclusion,  when  considering  estimation  with  delay,  it  must 
be  decided  if  the  improvement  in  the  estimate  is  worth  the  delay  in 
receiving  the  estimate  and  the  increased  complexity  of  the  optimum 
filter.  The  examples  in  the  following  section  will  serve  to  illustrate 
some  of  the  ideas  of  this  chapter. 

D.  EXAMPLES 

The  procedure  for  optimum  interp elation  with  delay  developed  in 
this  chapter  is  intimately  related  to  the  technique  for  optimum  filtering 
and  prediction  derived  in  Chapter  III.  In  order  to  clarify  this  rela¬ 
tionship,  the  examples  in  this  section  will  be  concerned  with  the  same 
two  stationary  stochastic  processes  that  were  used  in  Section  III-F.  to 
illustrate  the  ideas  of  Chapter  III.  For  the  first-order  Markov  process, 
the  optimum  interpolation  will  be  determined  for  a  delay  of  one  sample 
period  and  two  sample  periods.  For  the  stationary  process  with 


-  70  - 


exponential  cosine  autocorrelation,  the  optimum  interpolation  will  "be 
determined  for  a  delay  of  one  sample  period.  With  a  delay  of  a  whole 
number  of  sample  periods  &  is  equal  to  zero,  and  from  Eqs.  (l88)  and 
(202), 


Vv*) 


a2(V*> 


=  0 


6=0 


6=0 


(209) 


1.  Exponential  Autocorrelation 

The  stationary  first-order  Markov  process  has  the  autocorrelation 

3:xx(t)  =  e"P|j|  <21°) 

The  subscript  1  will  not  be  used  on  the  variables  in  this  example,  so 
the  expected  value  of  the  transition  matrix  will  be  denoted  with 


V  =  E[exp  (-P  Tn)] 


(211) 


and  the  variance  of  the  multiplicative  noise  will  be  denoted  by  r. 

The  model  of  the  optimum  filter  for  interpolation  with  a  delay  of 
both  one  and  two  sample  periods  is  presented  in  Fig.  10.  That  part  ■ef 
the  filter  within  the  dashed  rectangle  duplicates  the  model  of  the  opti¬ 
mum  filter  for  prediction  presented  in  Fig.  6  and  discussed  in  Section 
III-F.  The  weighting  k(t^t^+1)  is  a  scalar  which  is  determined  from 
Eq.  (197)  to  be 


k^td,td+l^ 


3  pW 


(212) 


The  weighting  k(td>td+g)  determined  from  Eq.  (199)  as 
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OPTIMUM  ESTIMATE 


FIG.  10.  MODEL  OF  THE  OPTIMUM  FILTER  FOR  ESTIMATION  WITH  DELAY  FOR 
EXAMPLE  1. 


k(td’td+2) 


^2p(td/td)  k(Wi) 

jl!  )+r 


(213) 


d+2'  d+l' 


and  after  substituting  Eq.  (212)  for  k(td)td+1)  in  Eq.  (213)  this 
can  be  written 


k(td,td+2) 


^P(td/td) 


1  - 


P(td+l/td) 


(214) 


The  covariance  P(tn+1An)  and  PAnAn)  are  given  by  Eqs.  (144)  and 
(l45)  as 


p(t  A  )  = 
x  n+1'  n' 


f 


2  - 


p(tn/tn-l)+r 


+  i  -  jir 


-  r 

P  (tn  An  )  =  - 

p<Wi>« 


(215) 
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The  steady-state  value  of  the  covariance  of  the  error  with  a  delay  of 
one  sample  period  p(t  /t  )  is  determined  from  Eq.  (203)  after 
substituting  Eq.  (212)  for  k(t^/t^+^)  and  the  result  is 


p(t  /t  )  =  p(t  /t  ) 
n-1  n'  ■  n  n 


1 


02p(t  /t  ) 
r  v  n'  n 


(216) 


The  steady-state  value  of  the  covariance  of  the  error  with  a  delay  of 
two  sample  periods  p(tQ  g/t  )  is  determined  from  Eqs.  (207)  and  (2l6) 
after  substituting  Eq.  (l^3)  for  k(t^)  and  Eq.  (2l4)  for  k(td  /t^  ) 

and  the  result  is 


^2p(tn/tn)  ^p(tn/tn)2r2 

^n+A)+? 


(217) 

The  steady-state  value  of  the  three  covariances  p(tQ/tn), 
p(tn  1/tn)j  and  p(tn  2 An)  are  Pl°tte(i  in  F1-g-  H  as  a  function  of 
the  variance  of  the  multiplicative  noise  r  for  ■ 


0  =  0.9091, 


(218) 


which  is  the  expected  value  of  the  transition  matrix  when  the  sample 
period  has  an  exponential  distribution  (purely  random  sampling)  with 
the  parameters  given  by  Eq.  (219). 


H  =  1.0 

P  =  0.1.  (219) 


Notice  how  delaying  the  estimate  decreases  the  covariance  of  the 
error 

2.  Exponential  Cosine  Autocorrelation 

The  second  example  is  the  stationary  process  with  exponential 
cosine  autocorrelation 
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VARIANCE  Of  THE  NOISE  T 


FIG.  II.  COVARIANCE  OF  THE  ERROR  FOR  EXAMPLE  1. 


3 


X1X1 


(t) 


cos  yT 


P?  sin  7  | 
(2P2  +  72) 


(220) 


which  was  discussed  in  Sect.  II-C  and  III-F.  For  a  delay  of  one  sample 
period,  the  weighting  matrix  is  determined  from  Eq.  (197)  to  be 


K<Vld+i>  ■ 


^lAl^d^d^  +  ^12P12^td^td^ 

Pll^d+lV4* 

^llP21^td^td^  +  ^12P22^td/'td^ 
Pll(td+l/td)+7 


(220) 


The  elements  of  the  covariance  matrix  Ptt^/t^)  can  be  determined 
from  Eq.  (l27)  as 
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The  steady-state  covariance  of  the  error  p, , (t  ,/t  )  for  a  delay  of 

11  n-1  n 

one  sample  period  can  be  obtained  in  terms  of  the  elements  of  the 
covariance  matrix  P(tn+^/t  )  by  substituting  Eqs.  (220)  and  (221)  into 
Eq.  (203),  and  the  result  is 


Pll(tn-l/tn)  =  T 


p11(t  ,/t  )r 

*11'  n+1  n' 


r2 


Pll(tn+l/tn)+r  tPu(tn+AJ+rl2 


2  2 


j?,  .  .  p,  .(t  , /t  )p,  .(t  ,/t  )  \  .  (222) 

li  lj  li'  n+1'  n  lj  n+l'  r*  ' 


i=l  j=l 


The  steady-state  value  of  the  two  covariances  p, , (t  /t  )  and 

11  n'  n 

p^(t  ^/t  )  are  plotted. in  Fig,  12  as  a  function  of  mhe  multiplicative 
noise  when  the  sample  period  has  a  geometric  distribution  (purely  random 
sampling)  with  the  parameters  given  below,  in  Eq.  (223) 


p  =  1.0 

P  =  0.1 

7  =  1-0,  (223) 

which  means  that  the  values  of  the  matrices  are  given  by  Table  4  in 
Sect  III-F.  Delaying  the  estimate  for  one  sample  period  decreases 
the  covariance  of  the  error,  but  not  quite  as  much  as  in  the  previous 
example . 


VARIANCE  OF  THE  NOISE  r 


FIG.  12  COVARIANCE  OF  THE  ERROR  FOR  EXAMPLE  2. 


V.  STATIONARY  STOCHASTIC  PROCESSES  WITH  RANDOM  SAMPLE  PERIODS 


A.  INTRODUCTION 

The  purpose  of  this  chapter  is  to  show  how  the  techniques  developed 
in  this  investigation  can  be  applied  to  a  problem  which  has  been  con¬ 
sidered  in  the  literature.  A  stationary  stochastic  process  with  a 
rational  power  spectrum  is  sampled  so  that  the  sample  period  is  a 
stationary  random  variable  with  a  known  first-order  probability  distri¬ 
bution.  There  is  no  multiplicative  noise;  therefore,  the  only  random 
parameter  of  the  process  is  the  randomness  in  the  sample  period.  It  is 
assumed  that  the  sample  points  are  being  observed  in  real  time;  there¬ 
fore,  although  the  sample  period  T^  is  a  random  parameter,  the  value 
of  the  random  parameter  Tn  can  be  determined  exactly  as  soon  as  the 
sample  points  tQ  and  tn+^  are  observed. 

In  Sect.  V-B,  the  continuous  optimum  estimate  of  the  process  is 
presented  when  the  sample  period  can  be  considered  as  a  known  time- 
varying  parameter.  This  optimum  estimate  is  generated  by  a  linear 
dynamic  filter  with  a  time-varying  matrix- valued  gain.  If  the  optimum 
estimate  can  be  delayed  'until  an  on-line  computer  iteratively  determines 
the  gain  at  each  sample  point,  then  the  optimum  estimate  can  be 
implemented. 

When  it  is  not  practical  to  recompute  the  gain  at  each  sample 
point,  the  linear  filter  that  generates  the  optimum  estimate  is  modified 
so  that  the  matrix-valued  gain  is  a  constant.  The  constant  value  chosen 
is  that  value  which  minimizes  the  trace  of  the  covariance  matrix  of  the 
error  in  the  limit  as  the  number  of  observations  n  approaches 
infinity. 

In  the  approach  to  the  problem  that  has  been  considered  in  the 
literature  [Refs,  13,1^]  transform  techniques  are  used  to  design  the 
linear  time -invariant  filter  that  minimizes  the  mean  square  error  of  the 
continuous  estimate.  In  Sect,  V-C  this  time-invariant  filter  is  com¬ 
pared  with  the  optimum  filter  for  the  case  in  which  the  stochastic  pro¬ 
cess  is  a  stationary  first-order  Markov  process. 
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B.  THE  OPTIMUM  ESTIMATE 


When  the  sample  period  of  a  stationary  stochastic  process  is  a 
time-varying  parameter  that  is  known  exactly,  the  optimum  current 
estimate  is  given  by  Eq.  (119)  which  can  be  written 

S(tA>  -  A-i’A-A-AAiA)-'4  A-i’A-A-ih  (221,> 

where  the  transition  matrix  (r)  is  given  by  Eq.  (17)  as 

®  (t)  =  eFx.  (225) 

The  matrix  block  diagram  in  Fig.  13  represents  the  linear  filter  that 
implements  the  optimum  estimate  in  Eq.  (224).  That  part  of  the  filter 
within  the  dashed  rectangle  duplicates  the  model  of  the  process  and  has 
the  impulse  response  $  (t).  The  switch  on  the  left,  representing  the 
sampling  operation,  and  the  :  ynchronous  switch  in  the  feedback  loop  of 


OPTIMUM 

PREDICTION 


FIG.  13.  MODEL  OF  THE  OPTIMUM  FILTER. 


the  filter  close  simultaneously,  so  that  the  sum  of  their  outputs  is 
zero  except  at  a  sample  point  t  when  the  sum  is  an  impulse  with  area 

y(tn)  -  ^(T^mt^/t^).  (226) 

The  weighting  matrix  K(tn)  is  determined  from  Eq.  (120)>  which 
is 


K(t  )  =  P(t  /t  ,  )MT<  M  P(t  /t  )  M' 
n  v  n'  n-1  \  v  n'  n-1' 


(227) 


The  covariance  P(t  /t^  )  is  calculated  iteratively  from  the  nonlinear 

difference  equation  (125),  which  can  be  written 

p<Vl'V  ■ 

V-  T 

/ft 

0(Tn-CT)GVGT  0T(Tn-cr)da,  (228) 
o 


where  the  substitution  has  been  made  from  Eq.  (48)  that 

Tn 

Q(t  )  =  f  <t>(Tn-CT)GVGT  <tT(Tn-o-)do-.  (229) 

o 


From  Eqs.  (227)  and  (228)  the  covariance  of  the  error  of  the  optimum 
estimate  during  any  time  in  the  sample  period  is 


p(tn+T/tn)  =  ®(t)  (X  -  K(tn)M^  P(tn/tn_1)  1T(r) 


I 

+  l~  d’(T-cr)  GVGT  'J>T(t-o-)  do. 


(230) 


The  weighting  K(t^)  for  the  optimum  estimate  in  Eq.  (227)  can 
be  determined  in  real  time  by  an  on-line  computer  if  a  slight  delay  is 
permitued  in  the  estimate.  The  following  procedure  must  be  repeated 
at  each  sample  point: 
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1.  Determine  the  sample  period  Tn  as  soon  as  the  sample  point 
tn+i  is  observed. 

2.  Calculate  and  store  the  covariance  of  the  error  P(tn+q/tn)  from 
Eq.  (228)  using  the  known  sample  period  Tn  and  the  stored 
covariance  P(tn/tn.p). 

3.  Calculate  the  weighting  matrix  K(tn+p)  and  set  the  gain  of  the 
filter  accordingly. 

The  idea  of  using  an  on-line  computer  to  calculate  the  weighting 
matrix  at  each  sample  point  did  not  originate  from  this  investigation. 

In  particular,  it  follows  directly  from  the  work  of  Kalman  [Ref.  1] 
when  his  results  are  carefully  applied  to  this  problem.  The  advantage 
to  this  procedure  is  that  the  optimum  estimate  is  obtained,  but  the 
corresponding  disadvantage  is  that  an  on-line  computer  must  be  available, 
and  the  optimum  estimate  must  be  delayed  until  the  computer  has  com¬ 
pleted  the  necessary^  calculations  to  determine  the  weighting  matrix. 

Sometimes  it  is  impractical  to  recalculate  the  weighting  matrix 
at  each  sample  point.  It  would  be  desirable  to  find  some  constant 
value  for  the  matrix  that  would  minimize  the  average  value  of  the  trace 
of  the  covariance  matrix  of  the  error.  The  constant  value  of  the  matrix 
which  minimizes  the  average  value  of  the  trace  in  the  limit  as  the 
number  of  observations  n  approaches  infinity  can  be  determined  by  the 
techniques  developed  in  Chapter  III. 

Assume  that  the  sample  period  Tn  is  an  unknown  random  parameter, 
rather  than  a  known  time-varying  parameter.  The  series  of  weighting 
matrices  K(t^)  can  determined  to  minimize  the  average  value  of  the 
trace  of  the  covariance  matrix  of  the  error.  In  Sect.  III-E  it  was 
shown  that  when  all  the  statistics  of  the  processes  were  stationary,  then 
the  optimum  estimate  would  become  a  stationary  dynamic  system  in  the 
limit.  This  means  that  in  the  limit  the  weighting  matrices  K(tn) 
will  approach  the  desired  constant  value  K 

The  weighting  matrix  K(t^)  is  determined  from  Eq.  (227);  thus 

K(tr)  .  "t(M  5(Wl>  “T)  _1  (231> 

The  average  value  of  the  covariance  P(^nAn  is  calculated  iteratively 
from  the  average  value  of  Eq  (228)  which  can  be  written 
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X 


-  EtMTn)1  #tM  P(tn/tn.i)  “'r1  «] 

^VVl^®  T(Tn>1+E[$  (T^VVl^ 

(232) 

where  the  covariance  Q(tn)  is  given  l>y  Eg.  (229)  and  all  the  expecta¬ 
tions  in  Eq.  (232)  are  with  respect  to  the  sample  period  T^.  The 
limiting  value  of  the  covariance 

Lim  P(tn+1/tn)  =  Lim  Pft/t^)  =  P  (233) 

n-m  n-«o 

can  be  obtained  by  iterating  Eq.  (232)  until  a  steady-state  solution 
is  obtained,  or  it  can  be  obtained  directly  from  the  steady-state 
solution  of  Eq.  (232).  The  desired  constant  value  of  the  weighting 
matrix  K  is  obtained  by  substituting  the  limiting  value  of  the 
covariance  P  into  Eq.  (231). 

The  filter  in  Fig.  13  with  a  constant  gain  K  is  not  the  optimum 
filter,  but  it  is  the  best  filter  from  a  sub-optimum  class  of  filters. 

In  the  following  section  the  filter  in  Fig.  13  is  compared  with  another 
filter  from  a  sub-optimum  class  of  filters--the  linear  time -invariant 
filter. 

C.  COMPARISON  WITH  THE  TIME -INVARIANT  FILTER 

When  it  is  desired  to  design  the  time -invariant  filter  with  infinite 
memory  that  minimizes  the  mean  square  error  of  the  estimate,  transform 
techniques  can  be  used  to  get  a  solution  to  the  problem  [Refs.  13,  lU]. 

In  the  transform  approach,  the  spectral  density  of  the  power  spectrum  of 
the  sampled  stochastic  process  is  first  determined  from  a  complex  convo¬ 
lution  integral,  which  can  be  evaluated  by  the  method  of  residues  in 
some  cases.  Then,  the  synthesis  procedure  is  based  upon  the  standard 
Wiener  spectral  factorization  of  the  sampled  power  spectrum.  The 
resulting  filter  is  continuous  and  time-invariant  and  gives  the  minimum 
mean  square  error  for  all  time  of  any  linear  filter  that  is  time  invariant. 
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Even  with  a  constant  gain,  the  filter  discussed  in  the  previous 
section  and  shown  in  Fig.  13  is  not  time  invariant.  It  contains  a  syn¬ 
chronous  switch  that  closes  at  the  sample  point,  so  that  the  character¬ 
istics  of  the  filter  depend  upon  the  sample  period.  A  simple  example 
will  serve  to  show  the  difference  between  the  two  filters.  The  perfor¬ 
mance  of  these  two  filters  will  be  compared  with  that  of  a  pure  stretcher 
or  time  varying  hold. 

Consider  the  stationary  first-order  Markov  process  discussed  in 
Sect.  II-B  with  the  autocorrelation  function 


The  process  is  sampled,  the  sample  period  is  a  stationary  random 
variable,  and  it  is  desired  to  reconstruct  the  original  process  from 
the  observed  random  variables  in  order  to  minimize  the  mean-square  error. 

The  optimum  filter  for  this  process  is  presented  in  Fig.  Ik.  The 
gain  of  the  filter  is  unity  for  any  probability  distribution  of  the 
sample  period.  The  filter  is  not  time  invariant  because  of  the 
synchronous  switch.  The  output  of  the  filter  is 


SAMPLING 

OPERATION 


FIG.  U.  OPTIMUM  FILTER. 


x(t  +t)  =  e‘^Tx(t  )  t  <  t  +  T  <t  , ,  (23*0 

n  '  n  n  —  n  n+1 

so  that  at  each  sample  point  the  error  in  the  estimate  is  zero.  During 
the  sample  period.,  the  error  is  due  to  the  random  input  to  the  process. 
The  mean -squared  error  due  to  the  covariance  of  the  random  input  is 
obtained  from  the  one -dimensional  version  of  Eq.  (229),  which  is 

p(tn+T/tn)  =  J  e‘^T'a)  23e'P^T*a^dCT  =  1  -  e2^J  tn-tn+T<tn+l  (235^ 

0 

where  the  substitution  has  been  made  from  Eq.  (9)  that 

vu  =  2P  (236) 

From  Eq.  (235)  the  average  over  all  time  of  the  covariance  of  the  error 
(or  the  mean-squared  error)  is 


-23T 

1  -  E[e  n] 

=  1  - (237) 

2P  E[Tn] 

The  linear  time-invariant  filter  for  the  stationary  first-order 
Markov  process  is  presented  in  Fig.  15.  The  impulse  response  of  the 
filter  g(r)  that  minimizes  the  mean-squared  error  is  dependent  on  the 
probability  distribution  of  the  sample  period  so  that  a  different  g(-r) 
must  be  determined  for  each  distribution.  The  best  time-invariant 
filter  and  the  corresponding  mean-squared  error  have  been  determined  in 
closed  form  for  two  probability  distributions  of  the  sample  period— the 
geometric  and  exponential  distributions.  The  probability  laws  governing 
these  two  distributions  have  already  been  discussed  in  Sect-.  II-E,  and 
the  pertinent  information  concerning  these  laws  is  summarized  in  Table  1 
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FIG.  15.  TIME-INVARIANT  FILTER. 


For  the  geometric  distribution  of  the  sample  period  (periodic 
sampling  with  independent  misses)  the  mean-squared  error  of  the  time- 
invariant  filter  is  given  by  Eq.  (238) 


Mean-squared  error  =  1 


(l-o-)(e2pT-l) 

_  y 

2PT 


(238) 


where  cx  is 


|(l-e' 


2PT 


)(1  -  f)  +  \\j  (l-e'2PT)2(l  -  f-f* 


4  2? 
P 


(l-e^) 


(239) 

and  p,  q,  and  T  are  parameters  of  the  distribution. 


For  the  exponential  distribution  of  the  sample  period  (purely 
random  sampling)  the  mean-squared  error  of  the  time-invariant  filter  is 
given  by  Eq.  (240). 


Mean- squared  error  = 


-  p  + 


+  2Pp 


P 


(240) 


where  p  is  a  parameter  of  the  distribution. 

It  is  interesting  to  compare  the  optimum  filter  and  the  time 
invariant  filter  with  a  simple  time  varying  hold  The  output  of  the 
time  varying  hold  is 
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X 


S(tn+T)  =  x(tn) 


t  <  t  +  T  <  t 

n  —  n  n+1 


The  mean  squared  error  of  an  estimate  at  time  t  +t  would  "be 


2(1  -  a2i3T) 


and  the  average  over  all  time  of  the  mean  squared  error  is 

T' 


E 


2(l-e^T)dT 


E[Tn] 


-PT 

=  2  -  2  (1  -  E[e  n]  1  /  P  E[Tn] 


In  Table  6  the  mean  squared  error  of  these  three  filters  is  compared 
for  period  sampling,  periodic  sampling  with  independent  misses,  and 
purely  random  sampling  for  p  equal  to  one-tenth.  The  expected  value  of 
the  sample  period  is  normalized  to  unity  and  the  value  of  the  other  para 
meters  of  the  probability  distribution  are  listed  in  Table  6. 


TABLE  6.  MEAN- SQUARE  EBROB  FOB  OPTIMUM  AND  TIME-INVARIANT  FILTERS. 


Purely  Bandoa  or  Poiaaon  Sampling 

- 

Periodic  Sampling  with  Indep 

endent  Miaaea 

Periodic  Sampling 

Probability  Diatribution 
o  f  S aap 1 e  Period 

Con  at  an  t 

Geoaetri e 

Exponen  ti al 

Paranetera  of 

T-l 

T-l/  2 

M-l 

Probability  Diatribution 

p-q-l/2 

Mean  Saaple  period 

T-l 

T/p-1 

1/M-l 

Variance  of  Saaple  Period 

0 

1/2 

1 

With  /3-0.1 

Mean-Square  Error 
for  Optiaua  Filter 

0.  OS  37 

0. 1311 

0.  1667 

With  >3-0 .  1 

Mean-Square  Error 
for  Tiae-Inrari an t 

Filter 

0.  0937 

0. 2458 

0.3583 

With  £-0.1 

Mean-Square  Error 
for  Tiae  Varying  Hold 

0.0967 

0.  1400 

0. 1111 
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For  periodic  sampling  the  optimum  filter  and  the  time  invariant 
filter  are  identical,  but  as  the  variance  of  the  sample  period  increases, 
the  mean  squared  error  of  the  time  invariant  filter  increases  much  more 
rapidly  than  that  of  the  optimum  filter.  The  mean  ^uared  error  of  the 
simple  time-varying  hold  is  only  slightly  greater  than  that  of  the  opti¬ 
mum  filter  for  periodic  sampling  and  also  when  the  sample  period  is  a 
random  variable.  In  this  example,  the  time  invariant  filter  does  not 
approximate  the  true  optimum  very  well  for  large  variance  of  the  sample 
period,  while  the  simple  time-varying  hold  does  approximate  the  true 
optimum.  This  shows  that  sometimes  a  "non-optimum"  time-varying  filter 
is  better  than  the  best  time  invariant  filter. 
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VI.  CONCLUSIONS 

A.  SUMMARY 

In  this  investigation  the  general  solution  has  been  derived  for 
the  problem  of  the  optimum  linear  estimation  of  a  sampled  stochastic 
process  with  random  parameters  that  can  be  adequately  approximated  by 
the  model  presented  in  Chapter  II.  The  random  parameters  are  independent 
from  one  sample  point  to  the  next  with  known  mean  and  covariance.  The 
solution  can  be  implemented  by  a  linear  dynamic  system  with  a  matrix¬ 
valued  gain  (or  gains)  calculated  iteratively  for  each  sample  point. 

For  high-order  complex  systems  these  computations  are  most  easily  per¬ 
formed  by  a  digital  computer. 

B.  SUGGESTIONS  FOR  FUTURE  WORK 

The  ideas  presented  in  Chapter  IV  represent  the  first  thorough 
investigation  of  optimum  interpolation  with  delay  for  a  nonstationary 
sampled  stochastic  process.  As  yet,  no  completely  satisfactory  theory 
exists  for  the  optimum  linear  interpolation  with  delay  for  a  non¬ 
stationary  continuous  process  It  should  be  possible  to  extend  the 
technique  developed  in  Chapter  IV  to  the  interpolation  with  delay  of  the 
continuous  stochastic  process  with  white  noise  added  to  the  measurements. 
For  filtering  and  prediction,  this  problem  was  formulated  by  Kalman  and 
Bucy  [Ref.  18]: 

Another  possibility  is  concerned  with  the  area  of  adaptive  systems. 
The  optimum  estimate  in  this  investigation  is  implemented  by  a  linear 
dynamic  system  that  uses  the  expected  value  of  the  random  parameters. 

The  only  place  the  covariance  of  the  random  parameters  is  needed  is  in 
the  iterative  calculation  of  the  matrix -valued  gain.  If  only  the 
expected  value  of  the  random  parameters  were  known,  then  the  techniques 
of  adaptive  systems  could  be  used  to  adjust  the  matrix- valued  gain  to 
react  to  measured  variations  in  the  covariance  of  the  random  parameters. 


-  87  - 


REFERENCES 


1.  Kalman,  R.  E.,  "A  new  approach  to  linear  filtering  and 

prediction  problems".  Journal  of  Basic  Engineering, 

Trans .  ASME,  Series  D,  82,  I960. 

2.  Bendat,  J.  S  ,  Principles  and  Applications  of  Random  Noise  Theory, 

John  Wiley,  New  York,  1958. 

3«  Sherman,  S.,  "Non -mean -square  error  criteria",  IRE  Trans,  on 
Information  Theory,  IT-4,  1958. 

4.  Gunckel,  T.  L. ,  "Optimum  design  of  sampled-data  systems  with 

random  parameters",  Stanford  Electronics  Laboratories, 
Technics!!  Report  No.  2102-2,  April  24,  1961. 

5.  Wiener,  N.,  The  Extrapolation,  Interpolation,  and  Smoothing  of 

Stationary  Time  Series,  John  Wiley,  New  York,  19^9* 

6.  Franklin,  G.  F.,  "The  optimum  synthesis  of  sampled-data  systems", 

Doctoral  Dissertation,  Department  of  Electrical  Engineering, 
Columbia  University,  1955- 

7.  Amara,  R.  C.,  "The  linear  least  squares  synthesis  of  continuous 

and  sampled  data  multivariable  systems",  Stanford  Electronics 
Laboratories,  Technical  Report  No.  40,  July  28,  1958. 

8.  Rosenbloom,  A.,  J  Heilfron,  D.  L.  Trautman,  "Analysis  of  linear 

systems  with  randomly  varying  inputs  and  parameters",  IRE 
Convention  Record,  Part  4,  1955- 

9.  Bergen,  A.  R.,  "Stability  of  systems  with  randomly  time-varying 

parameters",  IRE  Trans,  on  Automatic  Control,  AC-5,  4,  Septem¬ 
ber  i960. 

10.  Kalman,  R.  E,,  "Analysis  and  synthesis  of  linear  systems 

operating  on  randomly  sampled  data",  Doctoral  Dissertation, 
Dept,  of  Electrical  Engineering,  Columbia  University, 

August  10,  1957- 

11.  Kalman,  R.  E  and  R.  W.  Koepcke,  "Optimal  synthesis  of  linear 

sampling  control  systems  using  generalized  performance 
indexes",  ASME  Trans.,  Paper  No.  58-IRD-6,  1958,  p.  1821. 

12.  Kalman,  R  E, ,  "On  the  general  theory  of  control  systems", 

Proc.  of  the  First  IFAC  Moscow  Congress,  i960 . 

13-  Bergen,  A.  R  ,  "The  synthesis  of  optimum  random  sampling  systems", 
Technical  Report  T-T/133,  Electronics  Research  Laboratories, 
Columbia  University,  1957* 


-  88  - 


REFERENCES  (Continued) 


14.  Hsieh,  H,  C.,  "On  the  optimum  synthesis  of  random  sampling 

multipole  filters  with  stationary  inputs",  Report  No-  6l-8, 
University  of  California,  Department  of  Engineering,  Los 
Angeles,  California. 


15«  Beutler,  F.  J,,  "A  generalization  of  Wiener  optimum  filtering 
and  prediction".  Doctoral  dissertation,  California 
Institute  of  Technology,  Pasadena,  California. 

16.  Beutler,  F.  J.,  "Prediction  and  filtering  for  fandom  parameter 
systems",  IRE  Trans,  on  Information  Theory,  IT-4, 

December  1958 . 

17-  Shaw,  L.  G.,  "Dual  mode  filtering",  Stanford  Electronics 

laboratories,  Technical  Report  No.  2106-2,  September  26,  i960. 

18.  Kalman,  R.  E.  and  R.  S.  Bucy,  "New  results  in  linear  filtering 

and  prediction  theory",  Journal  of  Basic  Engineering,  Trans. 
ASME,  Series  D,  83,  1961. 

19.  Davenport,  W.  B.  and  W.  L.  Root,  An  Introduction  to  the  Theory 

of  Random  Signals  and  Noise,  McGraw-Hill,  New  York,  1958. 

20.  Coddington,  E.  A.  and  N.  Levinson,  Theory  of  Ordinary 

Differential  Equations,  McGraw-Hill,  New  York,  1955- 

21.  Gantmacher,  F.  R. ,  The  Theory  of  Matrices,  Vol.  1,  Chelsea, 

New  York,  1959- 

22.  Friedman,  B  ,  Principles  and  Techniques  of  Applied  Mathematics, 

John  Wiley,  New  York,  195&- 

23.  Parzen,  E.,  Modern  Probability  Theory  and  its  Application,  John 

Wiley,  New  York,  i960. 

24.  Doob,  J.  L  ,  Stochastic  Processes,  John  Wiley,  New  York,  1956. 

25.  Feller,  W.,  An  Introduction  to  Probability  Theory  and  its 

Applications,  Vol,  1,  John  Wiley,  New  York,  1950. 

26.  Greville,  T.  N.  E,,  "The  pseudoinverse  of  a  rectangular  or 

singular  matrix  and  its  application  to  the  solution  of 
systems  of  linear  equations",  SIAM  Review,  2,  1,  January 
I960. 


-  89  - 


t 


BASIC  DISTRIBUTION  LIST 
December  1961 


2 van 3  Signal  Laboratory 
Belmar,  New  Jersey 
1  Attn:  Head,  Solid-State  Dev. 
Br.  Dr.  H.  Jacobs 

Commanding  General,  USASRDL 
5  Ft.  Monmouth,  New  Jersey 
1  Attn:  SIGRA/SL-SC 

Chief,  West  Coast  Office 
1  Pasadena  2,  California 

Commanding  Officer,  Frankford 
1  Philadelphia,  Pa. 


Director,  NOL 
1  White  Oak,  Maryland 

Commande  r ♦  NADC 
Johnsvllle,  Pa. 

1  Attn:  NADC  Library 

Dept,  of  the  Army 
Washington  25,  D.C. 

1  Attn:  Research  Support  Div. 

Chief,  P&O  Division 
Washington  25*  D.C. 

1  Attn:  SIGOP-5 


Ballistics  Research  Laboratory 
Aberdeen  Proving  Grd,  Md. 

2  Attn:  L.A.  Dclsasso 

Commanding  General 
Ft .  Huachuca,  Arizona 
1  Attn:  Technical  Library 

Chief  of  Naval  Research 
Washington  25,  D.C. 

Z  Attn:  Code  427 

1  Attn:  Code  463 

Commanding  Officer,  ONR  Br.Off 
1  San  Francisco  9,  California 

Chief  Scientist,  ONR  BR.  Offic 
1  Pasadena,  California 

Comm.  Officer,  ONR  Br.  Office 
1  Chicago  1,  Illinois 


Office  of  the  Chief  ofEngr. 

1  Washington  25,  D.C. 

Hqs . ,  US AF  ( AFRDR-NU-3) 
Pentagon,  Washington  25,  D.C. 

1  Attn:  Mr.  Harry  Mulkey 

Chief  of  Staff,  USAF 
Washington  25,  D.C. 

2  Attn:  AFDRT-ER 

Commanding  Officer 
Aeronautical  Systems  Division 
Wright- Patterson  AFB, Ohio 
1  Attn:  ASRNEB 
1  Attn:  ASRNRS-2 

Executive  Director 
AF  Office  of  Scientific  Res. 
Washington  25,  D.r. 

1  Attn:  Code  SRYA 


Comm.  Officer,  ONR  Br.  Office 
1  New  York  13,  New  York 

Off leer- ln-Charge,  ONR 
15  New  York,  New  York 


1 

1 

1 

1 

1 

1 

1 

1 


Director,  NR  Laboratory 
Washington  25,  D.C. 
Attn:  Code  5240 
Attn:  Code  7100 
Attn:  Code  5430 
Attn:  Code  5200 
Attn:  Code  5300 
Attn:  Code  5400 
Attn:  Code  5270 
Attn:  Code  5210  GA 


AF  Special  Weapons  Ctr. 
Kirtland  AFB,  New  Mexico 
2  Attn:  SWOI 

Director,  AUL 
Maxwell  AFB,  Alabama 
1  Attn:  CR  4582 

Commander 

Patriot  AFB,  Florida 

1  Attn:  AFMTC  Tech .Lib. , MU-135 

Commanding  General 
Rome,  New  York 

2  Attn:  RCRW 


Chief,  Bureau  of  Ships 
Washington  25,  D.C. 

3  Attn:  Code  810 
1  Attn:  Code  8l6 
1  Attn:  Code  820 
1  Attn:  Code  840 


Commanding  General 
Bedford,  Massachusetts 
1  Attn:  CRTOTT-2,  Electronics 

AF  Office  of  Scientific  Res. 
Washington  25,  D.C. 

1  Attn:  Tech.  Library 


Commanding  Officer 
Diamond  Ordnance  Fuze  Labs . 
Washington  25,  D.C. 

2  Attn:  ORDTL  430,  Dr.  R.  Young 

U.S.  Dept,  of  Commerce 

1  Boulder,  Colorado 

2  Attn:  Miss  V.  Lincoln 

Director,  Nat.  Security  Agcy 
Ft.  George  G.  Meade,  Md. 

1  Attn:  REMP-2 

California  Inst,  of  Tech. 
Pasadena,  California 
1  Attn:  Prof.  L.M.  Field 

California  Inst,  of  Tech 
Jet  Propulsion  Laboratory 
Pasadena,  California 
1  Attn:  Documents  Library 

University  of  California 
Berkeley  4,  California 
1  Attn:  Prof.  R.M.  Saunders, CHir 

U.  Of  California,  Radiation  L. 

Berkeley,  Ca3 ifornia 
1  Attn:  Dr.  R,  K.  Wakerling 

Univ.  of  California 
Los  Angeles  24,  California 
1  Attn:  C.T.  Leondes 

Georgia  Institute  of  Tech. 
Atlanta,-  Georgia 

1  Attn:  Mrs.  J.H.  Crosland,  Lib. 

Harvard  University 
Cambridge  38,  Massachusetts 
1  Attn:  Mrs.  E.  Farkas,  Lib. 

Tl’’  noia  Inst,  of  Tech. 

Chicago  16,  Illinois 
1  Attn:  Dr.  Paul  C.  Yuen 
1  Attn:  Dr.  George  I.  Cohn 

University  of  Illinois 
Urbana,  Illinois 
1  Attn:  Prof.  D.  Alpert 
1  Attn:  Paul  Coleman,  EE  Dept. 

1  Attn:  Library  Serials  Dept. 

Johns  Hopkins  University 
Baltimore  2,  Maryland 
1  Attn:  Librarian 

Johns  Hopkiiu.  University 
Silver  Spring,  Maryland 
1  Attn:  Document  Library 


Chief,  Bureau  of  Aeronautics 
Washington  25,  D.C. 

1  Attn:  Code  RAAV-33 

Chief,  Bureau  of  Naval  Weapons 
Washington  25,  D.C. 

1  Attn:  RREN 
1  Attn:  RAAV-44 

Chief  of  Naval  Operations 
Washington  25,  D.C. 

1  Attn:  Code  0p-94T 

Comm.  Officer  and  Director 
1  San  Diego  52,  California 

U.S.  Naval  Post  Graduate  Sch. 
Monterey,  California 
1  Attn:  Tech.  Rpts .  Librarian 

Commander,  U.S.N. Missile  Ctr. 
Pt.  Mugu,  California 
1  Attn:  5321 

U.S.  Naval  Weapons  Lab. 
Dahlgren,  Virginia 
1  Attn:  Tech.  Library 

Naval  Ordnance  Laboratory 
Corona,  California 
1  Attn:  Library 


AF  Systems  Command 
Washington  25,  D.C. 

1  Attn:  RDRLT, Cap t .Zimmerman 
1  Attn:  RDRESE-Mr.  Presnell 

Hqa .  ARDC 

Washington  25,  D.C. 

1  Attn:  P.DRS 
1  Attn:  RDRC 

Assistant  Sec.  of  Defense 
Washington  25,  D.C. 

1  Attn:  Tech.  Library 

Director  of  Defense 
1  Washington  25,  D.C. 

U.S.  Coast  Guard 
Washington  25,  D.C. 

1  Attn:  E‘5E 

Dept,  of  Defense 
Washington  25,  D.C. 

1  Attn:  Code  7H,  Fred  Miller 

Advisory  Grp  on  Electron  Tubes 
1  New  York  13,  New  York 

Adv.  Grp  on  Reliability  of 
Elec.  Equipment 
1  Washington  25,  D.C. 


Mass.  Inst,  of  Technology 
Cambridge  39,  Mass . 

1  Attn:  Research  Lab.  of  Elec. 

Mass.  Inst,  of  Tech. 

Lincoln  Lab 
1  Lexington  73,  Mass . 

Director,  Cooley  Elecs .  Lab. 

1  Ann  Arbor,  Mich. 

Univ.  of  Michigan 
Ann  Arbor,  Michigan 
1  Attn:  Tech.  Documents  Serv. 

University  of  New  Mexico 
Albuquerque,  New  Mexico 
1  Attn:  R.K.  Moore,  Chm., 

Univ.  of  Nevada 
Reno,  Nevada 

1  Attn:  Dr.  Robt.  Manhart 

Northwestern  University 
Evanston,  Illinois 
1  Attn:  Walter  S.  Toth 

Ohio  State  University 
Columbus  10,  Ohio 
1  Attn:  Prof.  E.M.  Boone 


Oregon  State  College 
Corvallis,  Oregon 
1  Attn:  H.J.  Oorthuya 

Purdue  University 
Lafayette,  Indiana 
1  Attn:  Library 

University  of  Rochester 
Rochester  20,  New  York 
1  Attn:  Dr.  Gerald  H.  Cohen 

Stanford  Research  Institute 
Menlo  Park,  California 
1  Attn:  Documents  Center 

Villanova  University 
Villanova,  Pa. 

1  Attn:  Thomas  C.  Gabriele 

Yale  University 
New  Haven,  Connecticut 
1  Sloane  Physics  Laboratory 
1  Dept,  of  Elec.  Engr. 

1  Dunham  Laboratory,  Engr.  Lib. 


LEL, . Inc. 

Copiague,  L.I.,  New 
1  Attn:  Mr.  Robt.  S.  M*utner 


Lockheed  Aircraft  Corp^t^ 
Sunnyvale,  California 
1  Attn:  a.W.  Price,  Suv,8y3tenF 
1  Attn:  Dr.  Harris,  Mgp / 


The  M&xson  Elecs.  CorK™,*-^ 
New  York  1,  New  York  ‘>oratlon 
1  Attn:  M.  Simpson 


MITRE  Corporation 
Bedford,  Mass. 

1  Attn:  Constance  Wiling  . 

Monsanto  Corporation 
St.  Louis  66,  Missouri 
1  Attn:  E.  Orban,  Inorggpip 
Development 


National  Company,  Inc 
Maiden  48,  Mass. 

1  Attn:  Technical  LibrQ:r^an 


University  of  Puerto  Rico 
Mayaguez,  Puerto  Rico 
1  Attn:  Dr.  Braullo  E\ieno 

University  of  Saskatchewan 
Saskatoon,  Canada 
1  Attn:  Prof .N.F. Moody 

Uppsala  University 
Uppaala,  Sweden 
1  Attn:  Dr.  P.A.  Tove 

Admiral -Corporation 
Chicago  47,  Illinois 
1  Attn:  N.  Robertson,  Lib . 

Airborne  Instruments  Lab . 

Deer  Pk,  L.I.,  New  York 
1  Attn:  J.  Dyer,  VP 

Bell  Telephone  Labs. 

Murray  Hill,  New  Jersey 
1  Attn:  Elecs.  Res.  Dept. 

Bomac  Laboratories 
Beverly,  Mass. 

1  Attn:  Wellesley  J.  Dodds 

Columbia  Radiation  Lab. 

1  New  York  27 ,  New  York 

Cook  Research  Labs. 

1  Morton  Grove,  Illinois 

Cornell  Aeronautical  Labs., Inc 
1  Buffalo  21,  New  York 

General  Electric  TWT  Product 
Section 

Palo  Alto,  California 
1  Attn:  Technical  Llbrsry 


Nortronica  Tech. .Info  Aeencv 
1  Hawthorne,  California*  * 

Polarad  Corporation 
Long  Island  City  1,  New  York 
1  Attn:  A.H.  Sonnes che^n 
Chief  Systems 

RCA  Laboratories 
Princeton,  New  Jersey 
1  Attn:  Dr.  W.M.  Webster 
1  Attn:  Harwick  Johns 

Rarao -Wooldridge  Corpopat-ion 
Los  Angeles  45,  Californla 
1  Attn:  Dorothy  K.  Goas# 

The  Rand  Corporation 
Santa  Monica,  California 
1  Attn:  Margaret  Andersen  Lib. 

Raytheon  Company 
Burlington,  Massachusetts 
1  Attn:  Librarian 

Raytheon  Manufacturiru  co# 
Waltham,  Massachusetts 
1  Attn:  Librarian 

Roger  White  Electron  Devi  ren 
1  Haskell,  New  Jersey 

STL  Technical  Library 
Space  Tech.  Labs.  Inc. 

1  Los  Angeles  45,  California 

Sperry  Rand  Corporatfon 
1  Gainesville,  Florida 

Commanding  Officer 
1  Mountain  View,  Califo^ia 


General  Electric  Company 

Auburn,  New  York 
1  Attn:  Mr.  Gerald  C.  Huth 

General  Electric  Company 

Schenectady,  New  York 
1  Attn:  G.E.  Feiker 
1  Attn:  E.D .  McArthur 

Gilfillan  Brothers 

Los  Angeles,  California 
1  Attn:  Engineering  Library 

Hughes  Aircraft  Company 
1  Culver  City,  California 

Huggins  Laboratories 
1  Sunnyvs 1 e ,  Cal If o  m 1 s 

ITT  Federal  Laboratories 

Nutley  10,  New  Jersey 
1  Attn:  A.W.  McEwan 


Technical  Research  Gi>oUp 
1  Syosset,  L .  I . ,  New  Yo^ 

Texas  Instruments  Inc. 

1  Dallas,  Texas 
1  Attn:  Dr.  W.  Adcock 
1  Attn:  M.E.  Chun 


Varian  Associates 
Palo  Alto,  California 
1  Attn:  Technical  Lib*ary 


«cooiii(5uuuoc  Jiiecfcru  Corn 
Baltimore  3*  Maryland  p 


1  Attn: 


G.  Ross  Kilgore,  Mgr 
Applied  Res.  Iiept. 
Baltimore  Laboratory 


